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Abstract: In hard collisions at a hadron collider the most appropriate description of the 
initial state depends on what is measured in the final state. Parton distribution functions 
(PDFs) evolved to the hard collision scale Q are appropriate for inclusive observables, but 
not for measurements with a specific number of hard jets, leptons, and photons. Here the 
incoming protons are probed and lose their identity to an incoming jet at a scale \ib <C Q, 
and the initial state is described by universal beam functions. We discuss the field-theoretic 
treatment of beam functions, and show that the beam function has the same RG evolution as 
the jet function to all orders in perturbation theory. In contrast to PDF evolution, the beam 
function evolution does not mix quarks and gluons and changes the virtuality of the colliding 
parton at fixed momentum fraction. At ^b, the incoming jet can be described perturbatively, 
and we give a detailed derivation of the one-loop matching of the quark beam function onto 
quark and gluon PDFs. We compute the associated NLO Wilson coefficients and explicitly 
verify the cancellation of IR singularities. As an application, we give an expression for the 
next-to-next-to-leading logarithmic order (NNLL) resummed Drell-Yan beam thrust cross 
section. 
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1. Introduction 



The primary goal of the experiments at the LHC and Tevatron is to search for the Higgs 
particle and physics beyond the Standard Model through collisions at the energy frontier. 
The fact that the short-distance processes of interest are interlaced with QCD interactions 
complicates the search. A schematic picture of a proton-proton collision is shown in Fig. |]. A 
quark or gluon is extracted from each proton (the red circles labeled /), and emits initial-state 
radiation (I) prior to the hard short-distance collision (at H). The hard collision produces 
strongly interacting partons which hadronize into collimated jets of hadrons (Ji,2,3)) as well 
as non-strongly interacting particles (represented in the figure by the £ + £~). Finally, all the 
strongly interacting particles, including the spectators in the proton, interact with soft low- 
momentum gluons and can exchange perpendicular momentum by virtual Glauber gluons 
(both indicated by the short orange lines labeled S). 

The theoretical description of the collision is dramatically simplified for inclusive mea- 
surements, such as pp — > X£ + £~ , where one does not restrict the hadronic final state X. In 
this case, the cross section can be factorized as da = iJincl®/®/; where each / denotes a par- 
ton distribution function (PDF), which gives the probability of extracting a parton from the 
proton, while all other components of the collision are gathered together in a perturbatively 
calculable function H[ nc \. However, inclusive measurements do not necessarily contain all the 
desired information. Experimentally, identifying a certain hard-interaction process requires 
distinguishing between events that have a specific number of hard jets, leptons, or photons 
separated from each other and from the beam directions. Such measurements introduce new 
low energy scales and perturbative series with large double logarithms. For these situations 
it is necessary in the theoretical description to distinguish more of the ingredients in Fig. |], 
such as X, Jj, and S. Monte Carlo programs provide a widely used method to model the 
ingredients in the full cross section, da = H <S> f ® f <S>I ®X ®\\_ i Ji® S, using notions 
from QCD factorization and properties of QCD in the soft and collinear limits. Monte Carlos 




Figure 1: Schematic picture of a proton-(anti)proton collision at the LHC or Tevatron. 



- 2 - 



have the virtue of providing a general tool for any observable, but have the disadvantage of 
making model-dependent assumptions to combine the ingredients and to calculate some of 
them. For specific observables a better approach is to use factorization theorems (when they 
are available), since this provides a rigorous method of defining and combining the various 
ingredients. 

Here we investigate so-called beam functions, B = X f, which describe the part of 
Fig. |l] associated with the initial state. They incorporate PDF effects as well as initial- 
state radiation via functions X that can be computed in perturbation theory ffl. Below we 
will describe a particular class of measurements, which correspond to pp — > L + jets with 
L a non-hadronic final state such as Z — > £ + £~ or h — > 77. For these measurements, a 
rigorous factorization theorem has been proven that involves beam functions [Q]. We start 
by describing the general physical picture associated with beam functions, which suggests 
that they will have a wider role in describing cross sections for events with any number of 
distinguished jets, e.g. pp — > W/Z + n jets. That is, the beam functions have a more universal 
nature than what has been proven explicitly so far for the 0-jet case. 



The initial-state physics described by beam functions is illustrated in Fig. |2(a)| , and is 
characterized by three distinct scales ^a <C fis *C hh- At a low hadronic scale ^a the 
incoming proton contains partons of type k whose distribution of momentum is described 
by PDFs, /&(£', /M.)- Here £' is the momentum fraction relative to the (massless) proton 
momentum. Evolving fi to higher scales sums up single logarithms with the standard DGLAP 
evolution |, |, § §, § , 



%^) = E/y /*(£>)■ (1.1) 



This changes the type k and momentum fraction £' of the partons, but constrains them to 
remain inside the proton. At a scale (Ib, the measurement of radiation in the final state probes 



the proton, breaking it apart as shown in Fig. 2(a) and identifying a parton j with momentum 



fraction £ according to /xb). Measurements which have this effect at hb <C \in are those 
that directly or indirectly constrain energetic radiation in the forward direction, for example, 
by distinguishing hadrons in a central jet from those in the forward directions. The radiation 
emitted by the parton j builds up an incoming jet described by the function lij, and together 
these two ingredients form the beam function, 

B t (t',X,^ B ) = Z Zj x J^j(t',J,^B)m,^B)- (1.2) 

The sum indicates that the parton i in the jet need not be the same as the parton j in the 
PDF. The emissions also change the momentum fraction from £ to x and push the parton i 
off-shell with spacelike (transverse) virtuality —if < 0. The evolution for \x > hb sums up the 
double-logarithmic series associated with the t-channel singularity as if —> 0. It changes the 
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Figure 2: (a) Physics described by the beam function. Starting at a low hadronic scale /xa the proton 
is described by a PDF /. At the scale the proton is probed by measuring radiation in the final 
state, identifying a parton j described by fj(£,,(J-B)- Above /Kb, the initial state becomes an incoming 
jet described by Iij(t, x/£, fi) for an off-shell parton i with spacelike virtuality —t, which enters the 
hard interaction at [Ih- (b) Schematic picture of the final state for isolated Drell-Yan. 



virtuality t' of the parton i, while leaving its identity and momentum fraction unchanged, 

d f 

li—Bi(t,x,n)= dt'jU-t-t'i^Biit',*^)- (1.3) 



This evolution stops at the hard scale where the off-shell parton i enters the hard partonic 
collision. For [i > hb the initial state is also sensitive to soft radiation as shown by the orange 



wider angle gluons in Fig. 2(a) . For cases where the beam function description suffices this 
soft radiation eikonalizes, and the corresponding soft Wilson line is one component of the soft 
function S that appears in the factorized cross section. 

In general, a beam function combines the PDF with a description of all energetic initial- 
state radiation that is collinear to the incoming proton direction up to t <C Q 2 ■ The parton's 
virtuality t effectively measures the transverse spread of the radiation around the beam axis. 
The specific type of beam function may depend on details of the measurements, much as 
how jet functions depend on the algorithm used to identify radiation in the jet |0, H f|. 
Our discussion here will focus on the most inclusive beam function, which probes t through 
the measurement of hadrons in the entire forward hemisphere corresponding to the proton's 
direction. The utility of beam functions is that for a class of cross sections they provide a 
universal description of initial-state radiation that does not need to be modeled or computed 
on a case by case basis. 

An example of a factorization theorem that involves beam functions is the "isolated Drell- 
Yan" process, pp — >■ Xl + l~ . Here, as depicted in Fig. 2(b), X is allowed to contain forward 



energetic radiation in jets about the beam axis, but only soft wide-angle radiation with no 
central jets. The presence of energetic forward radiation is an unavoidable consequence for 
processes involving generic parton momentum fractions x that are away from the threshold 



- 4 - 



limit x — >• 1. There are of course many ways one might enforce the events to have no central 
jets. In Ref. Q, a smooth central jet veto is implemented by constructing a simple inclusive 
observable, called "beam thrust" , defined as 

e Y B+(Y) + e- Y B+(Y) x a E cra B+(Y) + x b E cm B b + (Y) 
TB = Q = ? • (L4) 

Here, q 2 and Y are the total invariant mass and rapidity of the leptons, Q = y/q 2 , and 

Q Y Q —Y /i r\ 

x a = -^e , x 6 = — — e (1.5) 

correspond to the partonic momentum fractions transferred to the leptons. The hadronic 
momenta Ba(Y) and B?(Y) measure the total momentum of all hadrons in the final state at 
rapidities y > Y and y < Y, respectively (where the momenta are measured in the hadronic 
center-of-mass frame of the collision and the rapidities are with respect to the beam axis). 
Their plus components are defined as B+(Y) = n a ■ B a (Y) and B^(Y) = n b ■ B b (Y) where 
ria = (1,0,0,1) and n b = (1,0,0,-1) are light-cone vectors corresponding to the directions 
of the incoming protons (with the beam axis taken along the z direction). The interpretation 
of beam thrust is analogous to thrust for e + e~ to jets, but with the thrust axis fixed to be 
the beam axis. For tb — 1 the hadronic final state contains hard radiation with momentum 
perpendicular to the beam axis of order Q, while tb <C 1 corresponds to two-jet like events 
with hard radiation of order Q only near the direction of the beams. The dependence of 
B^ b (Y) on Y accounts for asymmetric collisions where the partonic center-of-mass frame 
is boosted with respect to the hadronic center-of-mass frame. Requiring tb < ex.p(—2y B nt ) 
essentially vetoes hard radiation in a rapidity interval of size y^ — 1 around Y, i.e. in 
the region \y — Y\ < y B ut — 1, while radiation in the larger interval y B ut + 1 is essentially 
unconstrained, with a smooth transition in between. Thus, interesting values for y B ut are 
around 1 to 2. 

In Ref. |]l[|, a rigorous factorization theorem for the Drell-Yan beam thrust cross section 
for small tb was derived, 



da 



dq 2 dYdr B 



a Y^Hij{q 2 ,n) J dt a dt b Bi(t a ,x a ,Li) Bj(t b ,x b ,Li) 



-' (}Sb[Qtb - ta Q tb ^ 



, , Aqcd 



(1.6) 



using the formalism of the soft-collinear effective theory (SCET) [10, 11, O, [L3|], supplemented 



with arguments to rule out the presence of Glauber gluons partially based on Refs. [H, [l5| . 
The sum runs over partons ij = {uu,uu,dd, . . .}. The hard function Hij(q 2 ,fi) contains 
virtual radiation at the hard scale [iu — Q- It is given by the square of Wilson coefficients 
from matching the relevant QCD currents onto SCET currents (and hence it is identical to 
the hard function appearing in the threshold Drell-Yan factorization theorem). The beam 
functions Bi(t a , x a , yu) and Bj(t b , x b , fi) describe the formation of incoming jets prior to the 
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hard collision due to collinear radiation from the incoming partons, as described above. They 
are the initial-state analogs of the final-state jet functions Jj(i, fi) (appearing for example 
in the analogous factorization theorem for thrust in e + e~ —¥ 2 jets), which describe the 
formation of a jet from the outgoing partons produced in the hard interaction. In contrast 
to the jet functions, the beam functions depend on the parton's momentum fraction x in 
addition to its virtuality. For Eq. (|1.6|), the beam scale is set by \x% ~ ^/tbQ. Finally, the 
soft function Ss{k + } fi) describes the effect of soft radiation from the incoming partons on the 
measurement of tb, much like the soft function for thrust encodes the effects of soft radiation 
from the outgoing partons. It is defined in terms of incoming Wilson lines (instead of outgoing 
ones) and is sensitive to the soft scale fj,s — tbQ- 

The cross section for tb contains double logarithms In 2 tb which become large for small 
Tb < 1. The factorization theorem in Eq. ( [O] ) allows us to systematically resum these to 
all orders in perturbation theory. The logarithms of tb are split up into logarithms of the 
three scale ratios h/hh, h/hb, /VMS that are resummed by evaluating all functions at their 
natural scale, i.e. Hij at fin, Bi and Bj at fis, Sb at /j,s, and then RG evolving them to the 
(arbitrary) common scale ji. 

In this paper we give a detailed discussion of the beam function, including a derivation 
of results that were quoted in Ref. [IJ. We start in Sec. |2| with a discussion of several formal 
aspects of the quark and gluon beam functions, including their definition in terms of ma- 
trix elements of operators in SCET, their all-order renormalization properties, their analytic 
structure, and the operator product expansion in Eq. (|1.2[) relating the beam functions to 
PDFs. In particular, we prove that the beam functions obey the RGE in Eq. (|1.3[) with the 
same anomalous dimension as the jet function to all orders in perturbation theory. (A part 
of the proof is relegated to App. [B].) This result also implies that the anomalous dimension 
of the hemisphere soft functions with incoming Wilson lines is identical to the anomalous 
dimension of the hemisphere soft function with outgoing Wilson lines appearing in e + e~ — > 2 
jets. 

In Sec. ^, we perform the one-loop matching of the quark beam function onto quark and 
gluon PDFs. Using an offshellness regulator we first give explicit details of the calculations 
for the quark beam function and PDFs. We verify explicitly that the beam function contains 
the same IR singularities as the PDFs at one loop, and extract results for X qq and X qg at 
next-to- leading order (NLO). In App. we repeat the matching calculation for the quark 
beam function in pure dimensional regularization. 

Our results show that beam functions must be defined with zero-bin subtractions fl6[| , 
but that in the OPE the subtractions are frozen out into the Wilson coefficients Xjj. The 
subtractions are in fact necessary for the IR singularities in the beam functions and PDFs to 
agree. We briefly discuss why PDFs formulated with SCET collinear fields are identical with 
or without zero-bin subtractions. 

In Sec. ||, we first give the full expression for the resummed beam thrust cross section at 
small tb valid to any order in perturbation theory. The necessary ingredients for its evaluation 
at next-to-next-to-leading logarithmic (NNLL) order are collected in App. which are the 
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three-loop QCD cusp anomalous dimension, the two-loop standard anomalous dimensions, 
and the one-loop matching corrections for the various Wilson coefficients. (We also comment 
on the still missing ingredients required at N 3 LL order.) We then show plots of the quark 
beam function both at NLO in fixed-order perturbation theory and NNLL order in resummed 
perturbation theory. We also discuss the relative size of the quark and gluon contributions 
as well as the singular and nonsingular terms in the threshold limit x — ¥ 1. We conclude in 
Sec. § 

2. Beam Functions 
2.1 Definition 

In this subsection we discuss the definition of the quark and gluon beam functions in terms of 
matrix elements of operators in SCET, and compare them to the corresponding definition of 
the PDF. The operator language will be convenient to elucidate the renormalization structure 
and relation to jet functions in the following subsection. 

We first discuss some SCET ingredients that are relevant later on. We introduce light- 
cone vectors n 11 and with n 2 = n 2 = and n-n = 2 that are used to decompose four- vectors 
into light-cone coordinates p^ = (p + , p~ , p^_ ) , where p + = n ■ p, p~ = n ■ p and p^_ contains 
the components perpendicular to and n^. 

In SCET, the momentum p^ of energetic collinear particles moving close to the n direction 
is separated into large and small parts 

pt* =p» + p£ = n-p -j+Pn±+P?- ( 2 - X ) 

The large part p^ = (0,p~,pj_) has components p~ = n ■ p and p n ± ~ Ap _ , and the small 
residual piece pr = {pt iPV iPri) ~ p{^ 2 -> with A<C 1. The corresponding n-collinear 

quark and gluon fields are multipole expanded (with expansion parameter A). This means 
particles with different large components are described by separate quantum fields, £ n ,p(y) 
and A n ^p(y), which are distinguished by explicit momentum labels on the fields (in addition 
to the n label specifying the collinear direction). We use y to denote the position of the fields 
in the operators to reserve x for the parton momentum fractions. Two different types of 
collinear fields will be relevant for our discussion depending on whether or not they contain 

1/2 

perturbatively calculable components. For the beam functions A ~ , and the collinear 
modes have perturbative components with p 2 ~ Qtb 3> Aqcd- Collinear fields such as these 
are referred to as belonging to an SCETi theory. For collinear modes in the parton distribution 
functions A ~ Aqcd/Q and the collinear modes are nonperturbative with p 2 ~ Aq CD . These 
modes are a subset of the SCETj collinear modes and we will refer to their fields as belonging 
to SCETn. For much of our discussion the distinction between these two types of collinear 
modes is not important and we can just generically talk about collinear fields. When it is 
important we will refer explicitly to SCETj and SCETn. 
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Interactions between collinear fields cannot change the direction n but change the mo- 
mentum labels to satisfy label momentum conservation. Since the momentum labels are 
changed by interactions, it is convenient to use the short-hand notations 

Uv) = £ , A£(y) = £ A^(y) . (2.2) 

p^O p^O 

The sum over label momenta explicitly excludes the case p 11 = to avoid double-counting 
the soft degrees of freedom (described by separate soft quark and gluon fields). In practice 



when calculating matrix elements, this is implemented using zero-bin subtractions [16| or 



alternatively by dividing out matrix elements of Wilson lines [17, 18, 19]. The dependence on 
the label momentum is obtained using label momentum operators V n or Vl!± which return 
the sum of the minus or perpendicular label components of all n-collinear fields on which they 
act. 

The decomposition into label and residual momenta is not unique. Although the explicit 
dependence on the vectors and breaks Lorentz invariance, the theory must still be 
invariant under changes to and which preserve the power counting of the different mo- 
mentum components and the defining relations n 2 = n 2 = 0, n-n = 2. This reparametrization 
invariance (RPI) [p0| , ^] can be divided into three types. RPI-I and RPI-II transformations 
correspond to rotations of n and n. We will mainly use RPI-III under which and n^ 1 
transform as 

-> e a n" , n» -> e~ a n^ , (2.3) 

which implies that the vector components transform as p + — > e a p + and p~ — > e~ a p~ . In this 
way, the vector p^ stays invariant and Lorentz symmetry is restored within a cone about the 
direction of n^ 1 . Since Eq. ( [2 .3D only acts in the n-collinear sector, it is not equivalent to a 
spacetime boost of the whole physical system. 
We now define the following bare operators 

a g bare (y-,^) = e -^-/ 2 x n (y-^)|[^ -V n ) X n(0)] , 

6f™{y-,u) = e-^/ 2 tr{| Xn (y-^) [6(u - Pn)Xn(O)] } , 

O h g ^(y-,co) = -coe-^y-/ 2 B c n± ^(y-~) [5(u - V n )B^(0)] . (2.4) 

Their renormalization will be discussed in the next subsection. The corresponding renormal- 
ized operators are denoted as Oi(y~ , oj , fi) and are defined in Eq. ( [2.21 ) below. Here, p + is 



the momentum operator of the residual plus momentum and acts on everything to its right. 
The overall phase is included such that the Fourier-conjugate variable to y~ corresponds to 
the plus momentum of the initial-state radiation, see Eq. ( p.Sj ) below. The operator 5(u} — V n ) 
only acts inside the square brackets and forces the total sum of the minus labels of all fields 
in Xn(0) and B n j_(0) to be equal to to. The color indices of the quark fields are suppressed 
and summed over, c is an adjoint color index that is summed over, and the trace in Oq is 
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over spin. The operators are RPI-III invariant, because the transformation of the 5(lu — V n ) 
is compensated by that of the ijL in O q ^ and the overall to in O g . 
The fields 

Xn(y) = Wl{y) UV) , B n± = - g i W n(v) iD nJ^n{y)\ , (2-5) 

with i-D^j_ = + gA^ ± , are composite SCET fields of n-collinear quarks and gluons. In 
Eq. (|2.4|) they are at the positions y^ = y^n^/2 and y^ 1 = 0. The Wilson lines 



W n (y) 



E Q ^>{-=-n-A n (y) 

' n 



perms 



(2.6) 



are required to make Xn(y) and B^ ± (y) gauge invariant with respect to collinear gauge trans- 
formations [|ll], p^| . They are Wilson lines in label momentum space consisting of n-A n (y) 
collinear gluon fields. They sum up arbitrary emissions of n-collinear gluons from an n- 
collinear quark or gluon, which are 0(1) in the SCET power counting. Since W n (y) is lo- 
calized with respect to the residual position y, Xn(y) and B^ ± (y) are local operators for soft 
interactions. In SCETj the fields in Eqs. ( |2.4| ) and ( |2.5[) are those after the field redefini- 
tion decoupling soft gluons from collinear particles. Thus at leading order in the power 
counting these collinear fields do not interact with soft gluons through their Lagrangian and 
no longer transform under soft gauge transformations. Hence, the operators in Eq. ( |2.4[ ) are 
gauge invariant under both soft and collinear gauge transformations. The soft interactions 
with collinear particles are factorized into a soft function, which is a matrix element of soft 
Wilson lines. 

Note that our collinear fields in Eq. ([TJ) have continuous labels and hence are not the 
standard SCET fields with discrete labels. They only depend on the minus coordinate, y~, 
corresponding to the residual plus momentum, and not a full four- vector y^. As discussed 
in detail in the derivation of the factorization theorem in Ref . Q , it is convenient to absorb 
the residual minus and perpendicular components into the label momenta which then become 
continuous variables. For example, for the minus momentum (suppressing the perpendicular 
dependence) 

E e- ip ~ y+/2 xn,p-(y-,y + ) = E/ &PrZ-' l{r+p r )y+/2 Xn,Ay-) 

p- p- 

Ve- ip ^ +/2 Xn, P -(y")- (2.7) 



In this case, W n (y~n/2) can also be written in position space where all gluon fields sit 
at the same residual minus coordinate, y~ , and are path ordered in the plus coordinate 
(corresponding to the label minus momentum) from y + to infinity. 
Next, we introduce the Fourier-transformed operators 

0^ e (\cj\b + ,cj) = ±- [^e^y-^O^iy-^), (2.8) 

Z7T / Z\L0\ 
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and the corresponding renormalized operators Oi(\u\b + , u, /i) [see Eq. ( |2.22|) below]. For 
example, for the quark operator 

o h q -(\u\b\u) = i- /^^^(e^^ 

= Xn(0) 5(oob+ - LOp + ) | [6{u - Pn)Xn(0)] ■ (2.9) 

In the first step we used residual momentum conservation to shift the position of the field. 
Here we see that the overall phase in Eq. fl2,4| ) allows us to write the b + dependence in terms 
of 5{ub + — ujp + ), which means that b + measures the plus momentum of any intermediate 
state that is inserted between the fields. 

We divide by |u>| in Eq. ( [2.8D to make the integration measure of the Fourier transform 
RPI-III invariant. Using the absolute value \u\ ensures that the definition of the Fourier 
transform does not depend on the sign of ui and that the first argument of O q , t = \uj\b + , 
always has the same sign as b + . The Fourier-transformed operators are still RPI-III invariant 
and only depend on b + through the RPI-III invariant combination t. The beam functions are 
defined as the proton matrix elements of the renormalized operators Oi(t,uj,fi) in SCETi, 

B i (t,x = oj/p-,ii) = {p n {p-)\e{u)Oi{t,uj^)\p n (p-)) . (2.10) 

The matrix elements are always averaged over proton spins, which we suppress in our notation. 
Note that part of the definition in Eq. ( [2.10 ) is the choice of the direction n such that the 



proton states have no perpendicular momentum, P^ = P _ n^/2, which is why we denote them 
as \p n (P~)). By RPI-III invariance, the beam functions can then only depend on the RPI-III 
invariant variables t = cjb + and x = uj/P~. The restriction 6{ui) on the right-hand side of 
Eq. ( 2. 10| ) is included to enforce that the Xn(0), Xn(0), or B n ±(0) fields annihilate a quark, 



antiquark, or gluon out of the proton, as we discuss further at the beginning of Sec. |2.5| . 

The definition of the beam functions can be compared with that of the quark and gluon 
PDFs. In SCET, the PDFs are defined [Q] in terms of the RPI-III invariant operators 

Q^V) = <V)Xn(0)f [5(0/ -V n )Xn(0)] , 

Q^V) = 9(c/)tr{£ X n(0)[6{u/-V n )M0)]}, 

Q^ r V) = -^(wO K±,(0) [6(u/ - V n )B^(0)] , (2.11) 

as the proton matrix elements in SCETn of the corresponding renormalized operators Qi(u/, fi) 
defined in Eq. Q2.14Q below, 

fi{u//p-,n) = {p n {P-)\Qi{J,ii)\p n {P-)) . (2.12) 

By RPI-III invariance, the PDFs can only depend on the momentum fraction £ = u'/P~. 
Beyond tree level £ or to' are not the same as x or oj, which is why we denote them differently. 
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Without the additional 9{oj') in the operators in Eq. ( 2.11 ) the quark and anti-quark PDFs 
would combine into one function, with the quark PDF corresponding to oj > and the 
antiquark PDF to oj < 0. We explicitly separate these pieces to keep analogous definitions 
for the PDFs and beam functions. 



It is important to note that the SCETn collinear fields in Eq. ( [2.11 ) do not require zero 



bin subtractions, because as is well-known, the soft region does not contribute to the PDFs. 
If one makes the field redefinitions £ n — > Y£ n and A n — > YA n Y^ to decouple soft gluons, 
then the soft Wilson lines Y cancel in Eq. ( |2.11| ). Equivalently, if the fields in Eq. ( 2.11| ) 



include zero-bin subtractions then the subtractions will cancel in the sum of all diagrams, 
just like the soft gluons. (This is easy to see by formulating the zero-bin subtraction as a 
field redefinition [|l^] analogous to the soft one but with Wilson lines in a different light-cone 
direction.) In contrast, the SCETi collinear fields in the beam function operator in Eq. ( |2.4D 
must include zero-bin subtractions. We will see this explicitly at one loop in our PDF and 
beam function calculations in Sec. |3[ 

The SCET definitions of the PDFs are equivalent to the standard definition in terms of 
full QCD quark fields ip in position space. For example, the quark PDF in QCD is defined 
as 



m 



/ 9 (u//P- = 9(oj') j^^' y+l2 (pn(P-)\ [^(y + l)w(y+^)^m]y{P-)) ■ 

(2.13) 

The square brackets denote the renormalized operator. Here, the fields are separated along 
the n direction and the lightlike Wilson line W(y + n/2,0) is required to render the product 
of the fields gauge invariant. The relation to the SCET definition is that the SCETn fields 
m Eq. (pill ) (without zero-bin subtractions) involve a Fourier transform of ip in y + to give 



the conjugate variable oj'. The corresponding Wilson lines in Eq. ( |2.12 ) are precisely the 
W n contained in the definitions of Xn and S^ ± . Hence, the QCD and SCET definitions of 
the PDF are equivalent (provided of course that one uses the same renormalization scheme, 
which we do). 

Comparing Eq. ( |2 .11 ) to Eq. (2^), the difference between the beam functions and PDFs 



is that for the beam functions the fields are additionally separated along the n light-cone, 
with a large separation y~ S> y + corresponding to the small momentum b + <C oj. Thus, 
formulating equivalent definitions of the beam functions directly in QCD would be more 
challenging, as it would require QCD fields that are simultaneously separated in the n and 
n directions. For this case, it is not clear a priori how to obtain an unambiguous gauge- 
invariant definition, because Wilson lines connecting the fields along different paths are not 
equivalent. This ambiguity is resolved in SCETj, where the multipole expansion distinguishes 
the different scales and divides the possible gauge transformations into global, collinear, and 
soft transformations, allowing one to treat the separations along the two orthogonal light- 
cones independently. The large y~ separation corresponds to soft Wilson lines and soft gauge 
transformations that are independent from collinear gauge transformations corresponding to 
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the small y + dependence. As already mentioned, the operators in Eq. ( |2.4j ) are separately 
gauge invariant under both types of gauge transformations. 

2.2 Renormalization and RGE 

The beam functions and PDFs are defined as the matrix elements of renormalized operators. 
The renormalization of the operators immediately yields that of the functions defined by their 
matrix elements. In this subsection we derive the RG equations and show that the anomalous 
dimensions of the beam and jet functions are the same to all orders in perturbation theory. 

We start by considering the known renormalization of the PDF, but in the SCET operator 
language. The renormalized PDF operators are given in terms of the bare operators in 
Eq. ( [HI as 

3 

In general, operators with different i and u> can (and will) mix into each other, so the renormal- 
ization constant zf-ibj/oj' , fi) is a matrix in i,j and u, u'. RPI-III invariance then restricts the 
integration measure to be du'/uj' and zf-{uj /uj' , fi) to only depend on the ratio uj/uj' . Hence, 
the form of Eq. ( |2.14| ) is completely specified by the SCET symmetries. The [i independence 
of the bare operators Q^ are (cj) yields an RGE for the renormalized operators in MS 

j 

k 

where the inverse (Zf)~^{z } /j) is defined as 

Ej^( Zf )ik{^^) Z l j ^^) = ^S(l-z). (2.16) 
Taking the proton matrix element of Eq. ( |2.15| ) yields the RGE for the PDFs 

%/^)=E/f <{^)m'^)- ( 2 -!7) 

The solution of this RGE can be written in terms of an evolution function U* which acts on 
the initial PDF /j(£',jUo) an d takes it to /i(^, /u), 

m^) = (2-18) 

From Eq. ( ^17\) we can identify the anomalous dimensions 7^(2) in terms of the QCD 
splitting functions. For example, in dimensional regularization in the MS scheme, the one-loop 
anomalous dimensions for the quark PDF are the standard ones 

•>&(*, Ai) = ^^e(z)P qq (z), 7 / (z,n) = ^^e(z)P qg (z), (2.19) 
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with the q — )• qg and g — )• qq splitting functions 



9(1- 



z)(l + z 2 ) + ^(l 



^(l-zf + z 



9(1 -z) 



1 + z 2 



(2.20) 



The plus distribution Cq(x) = [9(x)/x] + is defined in the standard way, see Eq. ( A.2 ). For 
later convenience we do not include the overall color factors in the definitions in Eq. ( |2.20| ). 

We now go through an analogous discussion for the beam functions. The renormalized 
operators Oi(y~ , uj, /i) are given in terms of the bare operators in Eq. (|2.4|) by 



d bare (y - 



Z 



B 



,fj.)Oi(y ,u),fi) 



(2.21) 



where Z B (y~ /2\u\, [i) is the position-space renormalization constant. In App. |B|, we give 
an explicit proof that the beam function renormalization is multiplicative in this way to all 
orders in perturbation theory. 1 The underlying reason is that the renormalization of the 
theory should preserve locality, so renormalizing the nonlocal beam function operator should 
not affect the y~ separation between the fields. For example, mixing between operators 
with different y~ would destroy locality at distance scales within the validity range of the 
effective theory. RPI-III invariance then implies that Z B can only depend on the ratio y~ /2\oj\ 
(the factor of 1/2 is for convenience). In principle, one might think there could also be 



mixing between operators with different i or u in Eq. (2.21) [as was the case for the PDFs in 
Eq. ( p. 14 )]. Our derivation in App. ^ shows that this is not the case. 

Taking the Fourier transform of Eq. (2.21) according to Eq. (|2.8|), we find 



jdt' Z l B (t-t', fi) 



dy- 



(2.22) 



Since the bare operator is \x independent, taking the derivative with respect to /x, we find the 
RGE for the renormalized operator 



d 

fi-^-Oi(t,u,fi) 

IB (t,fl) 



dt'j B (t-t',»)Oi(t',u,»), 

jdt'(z B )-\t-t'^)^z B (t'^), 



(2.23) 



1 With our definitions of b + and t = they are always positive irrespective of the sign of w (i.e. for 

both beam and jet functions). Since y~ is the Fourier conjugate variable to b + , the Fourier-conjugate variable 
to t is u — y~ /2\ui\. The proof in the appendix, which is for oj > 0, shows that 7_b only depends on u through 
ln[i(w — i0)]. Its most general RPI-III invariant form is Jb(u, w/|w|). From Eq. (2.4) ^O^ arc (y~ , w)) = 
(0^ Te (— y~ , — u)) for any forward matrix element. Since the renormalization does not change the analytic 
structure, the same is true for the renormalized matrix elements so 7s(m, u}/\uj\) = 7b(— u, and also 

7s can only be a real function of w/|o;|. Because of its simple u dependence, we can conclude that 73 = 7b(m) 
only. With the tree-level boundary condition this then implies Zb = Zb(u). 
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where the inverse of Z B (t,n) is defined as usual, 

J at' {z B )-\t - 1', n) z l B (t', n) = s(t) . (2.24) 

Taking the proton matrix element of Eq. ( |2.23| ) we obtain the corresponding RGE for the 



beam function, 



IJL—Bi(t,x,^) = I dt'f B (t-t',n)Bi(t',x,n). 



(2.25) 



As discussed in App. [B], to all orders in perturbation theory the anomalous dimension has 
the form 



f B (t,^) = -2T l cnsp (a s ) -^Co 



+ Y B (a s )5(t) 



(2.26) 



where £q(x) = [0(x)/x]+ is defined in Eq. ( |A.2j ), (a s ) is the cusp anomalous dimension 
for quarks/antiquarks (i = q) or gluons (i = g), and j B (a s ) denotes the non-cusp part. Since 
there is no mixing between operators Oi(t, u, fj,) with different i or lj, the beam function RGE 
only changes the virtuality t but not the momentum fraction x and does not mix quark and 
gluon beam functions. By rescaling the plus distribution, 



Mo 



(2.27) 



1 C ( 1 

we can see that j l B (t,fj,) has logarithmic /i-dependence, which means that the RGE sums 
Sudakov double logarithms. 

The solution of the RGE in Eq. ( [2.25| ) with the form of the anomalous dimension in 
Eq. (p.26D is known |J, HI H- lt takes the form 



Bi(t,x,fi) = / dt' Bi(t - t',x,^ ) U B (t',Ho,n) 



where the evolution kernel can be written as 
U B (t,fi ,fi) 



Vb 



2 \ 2 



+ 6(t) 



(2.28) 



(2.29) 



T(l + rj B ) 

The distribution C v (x) is defined in Eq. (|A.2| ), and the RGE functions K B = K B ([iQ,[i) and 
rf B = r] B (fjLQ, fx) are given in Eq. ( D.7 ). 

The SCET quark, antiquark, and gluon jet functions are given by [^, 28] 

J q (up + +uj 2 L ,(j,) 



(2vr) 2 fdy 



2 \u)\ 



jp+y-/2 



tr(0 



L f [s(u + v n )6 2 (u± + V n± )Xn(0)] 



Jq(ujp + + CJ 2 , n) 

_ (2vr) 2 fdy 



(2vr) 2 



JV?-1 



0J 



dy_ 

2|w| 



Jp+y-/2 



(2.30) 
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where the notation [. . .]„ again denotes the renormalized operators. Here, we used the same 
conventions as for the beam functions where the large label momenta u and uj±_ are continuous, 
so the only position dependence of the fields is in the minus component. RPI invariance 
requires that the jet function only depends on the total invariant mass of the jet, p 2 = 
ujp + + When the jet function appears in a factorization theorem, the direction of the jet 
is either measured (e.g. by measuring the thrust axis in e + e~ — > 2 jets) or fixed by kinematics 
(e.g. in B — > X s 7 the jet direction is fixed by the direction of the photon) and n is chosen 
along the jet direction, so one typically has wj_ = 0. Taking the vacuum matrix element of 
0^ re (t,u>), we get 

^£(o\O h q ^(t,-u>)\o) 



)Xn(0) 



d 2 £± J„- are (t - u\) = J^ arc (t) = J^ are (t) . (2.31) 



In the last step we used that the quark and antiquark jet functions are the same. The 
analogous relation holds for the antiquark operator, Oq(t,UJ, fi). The u)± integral is bounded 
and does not lead to new UV divergences, because the jet function only has support for 
nonnegative argument, < uj± < t, and t is fixed. Similarly, for the gluon operator we have 

J^l-(0\O h g ^(t, -u,)|0) = J d 2 u ± jj-»(t - & ± ,p) = j^ c (t) . (2.32) 

The renormalization of J bare (t) does not depend on the choice of co± in Eq. ( |2.30| ). Since 
jW e (i) j g s i m piy an average over different choices for uj± it has the same renormalization. 
Hence Ji(t, jj) and Jj(t, //) have the same anomalous dimension, 

H-^Ji{t,fj,) = J d 2 uj±ds~/j(t -uj±- s,fi) Ji(s,fi) = J &t'ij(t - t',fjb) J d 2 u ± Ji(t' - 

, 'dt , 1 i j(t-t , , l i)J i (t , ,fx). (2.33) 



On the other hand, taking the vacuum matrix element of Eq. ( 2.23 ) we get 



d f 

^ Ji(t,fJ,) = J dt' f B (t - t' , fi) Ji(t' , fi) . (2.34) 

We thus conclude that the beam and jet function anomalous dimensions are identical to all 
orders in perturbation theory, 

7M*,M)=7^M)- (2-35) 
For the cusp part this result already follows from our explicit one-loop calculation, since 
is universal and its coefficients are the same at one loop. Our one-loop result provides a cross 
check for the identity of the one-loop non-cusp part of the anomalous dimension, which agree. 
Furthermore, jj(a s ) and hence ^ q B {a s ) can be obtained to three loops from Refs. [29, 3C], 
and for completeness the result is given in App. [p]. 
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2.3 Operator Product Expansion 

The difference between the beam function operators in Eq. ( [2.4D and the PDF operators 



in Eq. ( 2.11 ) is the additional separation in the y~ coordinate between the fields. Hence, 
by performing an operator product expansion about the limit y~ — > we can expand the 
renormalized operators Oi(y~ ,oj, fj) in terms of a sum over Qi(uj',fi), 



Oi(y ,u,fJ,) = Ji 



2w 



M 1 + 



3 



du/~ 



2|£J| 



(2.36) 



For completeness we included the identity operator on the right-hand side. The form of the 
matching coefficients 1^ and J{ is again constrained by RPI-III invariance so the structure 
of the OPE is completely determined by the SCET symmetries. Equation ( 2,36j ) encodes a 
matching computation between the operator Oi in SCETi, and the operators 1 and Qj in 
SCETn, where Ji and lij are the corresponding Wilson coefficients. 

Fourier transforming both sides of Eq. ( |2.36| ) with respect to y~ we get 



du/ 



Oi{t,uj,n) = Ji(t,iJ.)l + J2 I ^^'(t,^,^)Q J (w / , / u) + 0(^- 

* — * I CO \ UJ / \ U) 



y_ 



(2.37) 



Taking the vacuum matrix element of both sides, and using (0|Qj|0) = 0, we just get the 
coefficient of the identity operator on the right-hand side, which from Eqs. ( |2.31| ) and ( |2.32| ) 



is thus given by J{(t, fj,). Taking the proton matrix element of Eq. ( 2.37 ) with u > according 
to Eq. ( 2.10[ ), this first term drops out, because the jet functions only have support for — u > 
(or alternatively because the corresponding diagrams are disconnected), and we obtain the 
OPE for the beam function 



Bi(t,x,n) = X) / 



^1- 



l + 0(^) 



(2.38) 



For B g this equation was first derived in Ref. using a moment-space OPE for the matrix 
element (modulo missing the mixing contribution from the quark PDF). The higher-order 
power corrections in Eq. ( |2.38| ) must scale like 1/t and are therefore of 0(AQ CD /t) where 
Aq CD is the typical invariant mass of the partons in the proton. Equivalently, one can think 
of the scaling as (Aq CD /w)/6 + where Aq CD /u; is the typical plus momentum of the parton 
in the proton. These power corrections are given in terms of higher-twist proton structure 
functions. Since Eq. ( 2.37 ) is valid for t 3> Aq CD , this also means that we can calculate 
the matching coefficients in perturbation theory at the beam scale [i 2 B ~ t. This SCETj to 
SCETn matching calculation is carried out in the usual way by computing convenient matrix 
elements of the operators on both sides of Eq. ( |2.37j ) and extracting the Wilson coefficients 
from the difference. This is carried out at tree level in the next subsection, while the full 
one-loop matching calculation for the quark beam function is given in Sec. ||. On the other 
hand, for t ~ Aqcd the beam functions are nonperturbative and the OPE would require 
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an infinite set of higher-twist proton structure functions. In this case, the beam functions 
essentially become nonperturbative 6 + -dependent PDFs. 



The physical interpretation of the beam function OPE in Eq. ( |2.38| ) leads exactly to the 



physical picture shown in Fig. 2(a) and discussed in the introduction. At the beam scale 
Hb — t, the PDFs are evaluated and a parton j with momentum fraction £ is taken out 
of the proton. It then undergoes further collinear interactions, which are described by the 
perturbative Wilson coefficients z, //). By emitting collinear radiation it looses some of 
its momentum, and the final momentum fraction going into the hard interaction is x < £. 
In addition, the sum on j indicates that there is a mixing effect from terms without large 
logarithms, e.g. the quark beam function gets contributions from the quark, gluon, and 
antiquark PDFs. For example, when an incoming gluon from the proton pair-produces, 
with the quark participating in the hard interaction and the antiquark going into the beam 
remnant, then this is a mixing of the gluon PDF into the quark beam function. These are 
the physical effects that would usually be described by the PDF evolution. The difference is 
that once we are above the beam scale these effects only cause non-logarithmic perturbative 
corrections, which means the parton mixing and ^-reshuffling now appears in the matching, 
while the RG evolution of the beam function only changes t, as we saw above. In Sec. ||, we 
will see that these matching corrections are still important numerically and must be taken 
into account. For example, since the gluon PDF at small £ is very large compared to the 
quark and antiquark PDFs, it still gives an important contribution to the quark and antiquark 
beam functions. 

The consistency of the RGE requires that the fi dependence of the Wilson coefficients 
Iij(t,z,n) turns the RG running of the PDFs into the proper RG running of the beam 
functions. Taking the \i derivative of Eq. ( |2.38| ) we find the evolution equation for the Wilson 
coefficients 

H^-Iifaz,!*) =J2j dt '^r T ik(t-t',^,v) [i B (t\ii)8 kj 8{l -z') -Sit'^iz'^) 

k 

(2.39) 

The solution to this RGE can be easily obtained in terms of the evolution factors for the PDF 
and beam function in Eqs. (JH|) and ( jEgp , 



k 

Hence as expected, the RGE running of Xij(t, z, fi) cancels the running of the PDFs and adds 
in the running of the beam function. 

2.4 Tree-level Matching onto PDFs 

To illustrate the application of the OPE, we will calculate the Wilson coefficients Zjj at 
tree level, starting with X qq . We can use any external states for the computation of the 
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Figure 3: Tree-level diagram for the quark PDF (a) and the quark beam function (b). For the latter, 
the y~~ coordinate separation in the operator is indicated by drawing separated vertices for each field. 

Wilson coefficient as long as they have nonzero overlap with our operator. Thus, we pick 
the simplest choice, n-collinear quark and gluon states, \q n (p)} and \g n (p)}, with momentum 
= (p+,p~,0) where p~ > is the large momentum. In the following section we will use 
a small p + < as an IR regulator, but otherwise p + is set to zero. The tree-level diagrams 



with an external quark for the quark PDF and beam function are shown in Figs. 3(a) and 



3(b) , They give 



(qn(p)\Qq{u , p)\q n (p)) i0) = 9(uo')u n (p)5(J -p ) t;U n {p) = 9(u')5(l -Qj'/p ) , 

(q n (p)\O q (t,u>, v)\qn(p)) {0) = Un(p) S(t) S(oj - p~) K n {p) = 8(t) 6(1 - u>/p-) . (2.41) 

Here and in the following the superscript (i) indicates the 0(a l s ) contribution. Note that the 
results in Eq. ( p. 41 ) are the same whether we use a state with fixed spin and color or whether 



we average over spin and color. Taking the matrix element of both sides of Eq. ( 2.371 ) anci 
using Eq. ( 2.4l| ), we can read off the tree-level matching coefficient 



T<${t,z,y) =jg } (t,z,/x) = 5(t)5(l -z). (2.42) 
Similarly, the tree-level results for the gluon PDF and beam function are 

(gn(p)\Q g (u'^)\9n(p)) (0) = 9(U>')S(1 - J/p-) , 

(g n (p)\O g (t, oo, n)\g n {p)f ] = -oo s*-eS(t) 5(oo - p~) = 5(t) 5(1 - oo/p~) , (2.43) 
leading to 

l$(t,z, f i) = 6(t)6(l-z). (2.44) 

Finally, since at tree level the quark (gluon) matrix elements of the gluon (quark) operators 
vanish, 

(gn(p)\Q q (^,^\gn(p)f ] = {qn(p)\Q 9 (u'^)\q n (p)) m =0, 

(gn(p)\O q (t,L0,v)\gn(p)) {0) = (qn(p)\O g (t,co,fi)\q n (p)) {0) =0, (2.45) 
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we obtain 

lj%\t,z,»)=l$(t,z,») = 0. (2.46) 

To summarize, the complete tree-level results are 

l^(t,z,v) = 6^)6(1 -z), Bl 0) (t,x, t i) = 5(t)f i (x, f i). (2.47) 

The interpretation is simply that at tree level the parton taken out of the proton goes straight 
into the hard interaction. However, even at tree level the OPE already provides nontrivial 
information. From our general discussion we know that the matching should be performed 
at the beam scale fj, B ~ t to avoid large logarithms in the 0(a s ) terms, and this determines 
the scale at which the PDFs must be evaluated to be fi = hb- 

2.5 Analytic Structure and Time-Ordered Products 

In this subsection we discuss the analytic structure of the beam functions. For the OPE 
matching calculation we want to calculate partonic matrix elements of O q (t,to, fi). For this 
purpose it is convenient to relate the matrix elements of the products of fields in O q (t,oo,n) 
to discontinuities of matrix elements of time-ordered products of fields, since the latter are 
easily evaluated using standard Feynman rules. For notational simplicity we only consider 
the quark operator O q (t,oo) and suppress the spin indices and /x dependence. The discussion 
for the antiquark and gluon operators are analogous. 

We are interested in the forward matrix element of O q (t, oj) between some n-collinear state 
\Pn) = \Pn(p + ,P~)) with large momentum p~ and small residual momentum p + . Inserting a 
complete set of states ^2x\ X )(X\, we get 

( Pn \O q (t,u;)\p n )=J2(pn\M^p(t-HP + )\ X )( X \^ (2-48) 
X 

= Y J 5 ( t ~ \"\Px) S (" ~P~ +P X )(Pn\xn(0)^\x)(X\xn(0)\pn) ■ 
X 

The 5(oj — Vn) by definition only acts on the field inside the square bracket, returning its minus 
momentum, which by momentum conservation must be equal to the difference of the minus 
momenta of the external states. Since oj = p~ — p^, requiring to > implies p^ < p~ . This 
means that the action of the field reduces the momentum of the initial state so it effectively 
annihilates a parton in the initial state \p n ). Similarly, for to < we would have p% > p~ and 
the field would effectively create an antiquark in (X\. Also, since \X) are physical states, we 
have p^r > so oj < p~ and t = \oo\p% > 0. 

Hence, for the beam function, where \p n ) = \p n (P~)) is the proton state, the restriction 
to to > in its definition, Eq. ( [2.10 ), enforces that we indeed take a quark out of the proton. 



(Note that oj < does not correspond to the anti-quark beam function.) Taking the states 
\X) to be a complete set of physical intermediate states, the beam function has the physical 
support 

Q<x<1 _Pjmtn <1} t>oop+ min >0, (2.49) 
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where min > are the smallest possible momenta (which are strictly positive because with 
an incoming proton \X) can neither be massless nor the vacuum state). For the jet function 
the external state is the vacuum \p n ) = |0) yielding 5{uj + p x )> so the matrix element in 
Eq. ( [2.48 ) vanishes for uj > 0. 

Next, consider the following time-ordered analog of (p n \0(oj)O q (t, uj)\p n ^, 



(p n \T q (ujb + ',0j)\p n ) 



W) f^y_ e i(b+-p+)y-/2 {p 
2>tu 



2tt 



T{xn(y~^[S(u-Vn)Xn(0)]} 



Pn, 

(2.50) 



Writing out the time-ordering, 



T\ Xn[y 



2/2 

0(y-)Xn (jT|) | [S(U - V n )Xn(0)} ~ 0(-jT) [5{u> - V n )Xn(0)} Xnfy 



[6{u>-r n ) X n(o)]} 

_n\fi 
i — l — 



2/2 



using 



(2.51) 



(2.52) 



K + iO 

inserting a complete set of states, and translating the fields to spacetime position zero, we 
arrive at 



(p n \T q (ub + ,uj)\p n ) 



(27T) 



2uj K + iO 



E 

x 



X b+ - p +- K)y - / 2 6{cj _ p - +p - ) (p n \xn(0)^\x){X\ X n(0)\Pn) 



+ e i(b++p++K)y-/2 5 ( U +p -_ ( Pn | Xn(0 )|X>(x|xn(0)| 



Pn 



iO(u>) 
2-ituj 



E 

x 



5{uj-p +p x ) 
b+ -p x + iO 

_ 8(w +P~ ~P~x 
b+ + Px ~ i0 



Pr, 



Xn(0)^\x)(X\ Xn (0)\p n ) 



( Pn \Xn(0)\X)(x Xn(0)t 



Pn 



(2.53) 



The first term creates a cut in the complex b + plane for b + > p~x m i n - This cut is shown as 
the dark red line in Fig. ||. The second term produces a cut at b + < — p^- min , shown as the 
light blue line in Fig. ||. 

The beam function matrix element in Eq. ( |2.48| ) can be identified as precisely the dis- 
continuity of the first term in Eq. (2.53) with respect to b + . Thus, for the beam function we 
have 

B q (ub+,u) = Disc 6+>0 (p n (p-)\T q (ub + ,oj)\ Pn (p-)) . (2.54) 
Taking the discontinuity only for b + > ensures that we only pick out the cut due to the 



first term in Eq. (|2.53| ). Here, the discontinuity of a function g{x) for x > xq is defined as 

Disc^Tv, q{x) = lim 6(x - Xn) I q(x + iff) - a(x - iff) I , (2.55) 



. -x gix) = lim 6(x - xq) [g(x + iff) - g(x - iff)] , 
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Figure 4: Cuts in the complex b + plane for the time-ordered product in Eq. ( 2.53 ) 



and we used Eq. ( |A.7|) to take the discontinuity of l/(b + — pj) 

i 1 1 



Disc b + 



>o 



2tt \oj b + 



Px 



5(b + -p+) = 6(ub A 



(2.56) 



Since we explicitly specify how to take the discontinuity, we can drop the iO prescription in 
the denominators. (Alternatively, we could multiply by i and take the imaginary part using 
the iO prescription.) Since b + and t = \ui\b + always have the same sign we can also take the 
discontinuity for t > 0, so 



(p n \6(u))O q (t,u)\p n ) = Disc i>0 (p n \T q (t,u))\p n ) 



(2.57) 



For the matching calculation \p n ) is a partonic quark or gluon state. For any contribu- 
tions with real radiation in the intermediate state, i.e. diagrams where the two \n or B n ± 
fields in the operator Oi are joined by a series of propagators and vertices, we can use the 
standard Feynman rules to evaluate the time-ordered matrix element of T„(t, uS). However, 
with partonic external states, we can also have the vacuum state as an intermediate state, 
because the fields in the operator are spacetime separated. For such purely virtual contribu- 
tions it is simpler to directly start from O q (t,uj), insert the vacuum state between the fields, 
and then use standard Feynman rules to separately compute the two pieces (pn Xn (0)^/2 0) 
and (0|xn(0) Pn)- I n fact, this is exactly what we already did in our tree-level calculation in 
Sec. |2.4| , and we will see another example in Sec. ||. Thus, we will obtain the total partonic 
matrix element as 



(p n \e(uj)O q (t,u)\p n ) = ( Pn \8(u)O q (t,uj)\p n ) v . Ttual + {p n \6{uj)O q {t,Lj)\p n ) T ^ tion 

connected 

+ Disc t>0 <Pn|T g (t,^)|p n > connected . 



S(t)5(u-p-)( Pn x„(0)f 0) <0|xn(0)|p„> 

\ Z I connected 



connected 



(2.58) 



The virtual contribution must be kept, since it only looks superficially disconnected because 
the operator itself is spacetime separated. As always, we still disregard genuinely disconnected 
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diagrams, e.g. diagrams involving vacuum bubbles, when calculating the matrix elements in 
the second line. 



3. NLO Calculation of the Quark Beam Function 

In this section, we compute the matching coefficients I qq (t,z,fj,) and I qg (t,z,/j,) in the OPE 
for the quark beam function in Eq. ( [2.38 ) to next-to-leading order in a s (/x). As explained in 



Sec. and Sec. |2.4| , this can be done by computing the partonic matrix elements of both 
sides of Eq. Q2.37| ) to NLO. We use the same n-collinear quark and gluon states, \q n ) = \q n (p)) 



and \g n ) = \g n (p)), as i n the tree-level matching in Sec. 2.4 , with momentum p^ = {p + ,p~ , 0). 
Since only I qq (t, z, jj) is nonzero at leading order, we will only need the NLO matrix elements 
of the quark operators, O q (t,uj,n) and Q q (uj,fj,). We write the results for all matrix elements 
in terms of the RPI-III invariant variables (in this section we will always have uj > 0) 

t = ub + , t' = —ojp + = —zp + p~ , z = — . (3-1) 

p 

Here, z is the partonic momentum fraction of the quark annihilated by the operator relative 
to the momentum of the incoming quark or gluon, and will coincide with the argument of 
lij(t,z,n). 

To regulate the UV we use dimensional regularization with d = 4 — 2e dimensions and 
renormalize using the MS scheme. Since the matching coefficients in the OPE must be IR 
finite, the matrix elements of O q and Q q must have the same IR divergences, i.e., the beam 
function must contain the same IR divergences as the PDF. To explicitly check that this is 
the case, we separate the UV and IR divergences by regulating the IR with a small p + < 0. 
This forces the external states to have a small offshellness p + p~ < 0, and since p + p~ = —t'/z 
the IR divergences will appear as hit'. This also allows us to directly obtain the one-loop 
renormalization constants and anomalous dimensions for O q and Q q from their one-loop 
matrix elements. 

We first compute the renormalized one-loop matrix elements of the quark PDF operator 



Q q in Sec. 3.1. This calculation of the PDF for general x using the SCET operator definition 



and with an offshellness IR regulator is quite instructive, both by itself and in comparison to 



the beam function calculation, which is why we give it in some detail. In Sec. |3.2j , we compute 
the renormalized one-loop matrix elements of the quark beam function operator O q . Finally 
in Sec. 3^3, we use these results to extract expressions for I qq (t, z, fi) and I qg (t, z, fj) valid to 



NLO. 

Assuming that the IR divergences in the beam function and PDF will cancel, the matching 
calculation can be performed more easily using dimensional regularization for both UV and 
IR. We do this as an illustrative exercise in App. ^| which, as it should, yields the same result 
for the matching coefficients. 
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Figure 5: Nonzero one-loop diagrams for the quark PDF. The minus momentum u enters the vertex 
through its outgoing fermion line and leaves through its incoming fermion line. Diagram (c) represents 
the inclusion of the wave-function renormalization constant for the renormalized fields together with 
the corresponding residue factor in the LSZ formula for the S'-matrix. Diagrams (b) and (c) have 
symmetric counterparts which are included in their computation. 

3.1 PDF with Offshellness Infrared Regulator 

We start by calculating the bare S'-matrix elements 

(qn( P )\Q h q arc n\q n (p)) , {9n{p)\Q h r^)\g n {p)) , (3.2) 

using Feynman gauge to compute the gauge- invariant sum of all diagrams. The relevant one- 
loop diagrams are shown in Fig. [|. Since Q q is a local SCET operator, we can use the usual 



time-ordered Feynman rules in SCET (without any of the complications discussed in Sec. 2.5 
for O q ). The collinear q n q n g n vertex factor is 

ig I*V£(p, t) | with V*(p, I) = n» + t^L + Zpi _ t±t±^ , (3.3) 

where p^ and -P are the label momenta of the outgoing and incoming quark lines. (Because we 
have a single collinear direction the computation can also be done with QCD Feynman rules, 
still accounting for zero-bin subtractions, with the only difference being the Dirac algebra in 
the numerator of the loop integral. We checked that the final results for each diagram are 
indeed the same either way.) 
The diagram in Fig. |5(a)| is 

/ Inbare, n| \(«0 ./^V 2, f U n {p)V» {p, l)V n »(l, p) jujp) (T ) 2 

( qn \Q g H|g„) =-\—) 9CF jj^y {( 2 + iom _ p) 2 + iQ] ). 

_ (3-4) 
where g = g(fj-) is the renormalized MS coupling. The Dirac algebra for the numerator gives 

u n (p)V£(p,£)V nil (£,p)^u n (p)(r) 2 = u n (phlhhl±^n(p) = P~ (d - . (3.5) 

To compute the loop integral we write d d £ = d£ + d£ _ d d_2 ^j_/2, where £± is Euclidean, so 
£\ = — The £ + integral is done by contour integration as follows. For £~ < all poles are 
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above the axis and for t~ > p~ all poles are below the axis, so both cases give zero. Hence, 
the £~ integration range is restricted to < £~ < p~, where there is a double pole below 
the axis from the l/(£ 2 + iO) 2 and a single pole above the axis from the l/[(£ — p) 2 + iO]. 
Taking the single pole above amounts to replacing the second denominator by 2ni/(£~ — p~) 
and setting 



v 



(3.6) 



everywhere else. After performing the contour integral the iO have served their purpose and 
can be set to zero everywhere. The £~ integral is trivial using the 5(£~ — uS) and turns the £~ 
limits into an overall 9(u)9(p~ — uS). The remaining £± integration is done in d — 2 = 2(1 — e) 
Euclidean dimensions as usual. Putting everything together, we obtain 



(q n \Q 



bare / 



(U))\q t 

V 4vr 

a s (fi)C F 

2vr 
a s {p)C F 

2vr 



>(«) 

'g 2 C F e{u)9(p--u) 
6(z)9{l- z)T{e 



(d-2)(p--w) fd d ~ 2 i 



e^p 2 



4irp~ 

e 

(1- 



(27r) d - 2 [p + (1 - z ) t f 



(1 



9{z)9{l- z)(l- z) 



In- 



a 

J 2 



ln(l - z) - 2 



(3.7) 



where in the last line we expanded in e. 

In the diagram in Fig. 5(b), the gluon is annihilated by the Wilson line inside one of 
the Xn fields. The contraction with the one in Xn is oc 5(£~ — u) and the contraction with 
the one in Xn is oc S(p~ — uj). The 1/V n in the Wilson lines [see Eq. (|2.6| )1 contributes a 
factor \j((r — p~) with a relative minus sign between the two contractions. (There is also a 
diagram where the gluon connects both Wilson lines which vanishes because the Wilson lines 



only contain n-A gluons and we use Feynman gauge.) Adding Fig. 5(b) and its mirror graph, 
which gives an identical contribution, we get 



2i 



e^n 2 



Air 
a s (fi)C F 

7T 
7T 



g 2 c F 



d d £ 



(p)V n ^ Un (p)£ 



e^y 2 
—p + p~ 
e< E n 2 
t> 



(2ir) d (£' - p~){£ 2 + i0)[(£ - p) 2 + iO] 



d£- 9{£')9{p-- £~)-^ 
9(z)9(l - z)z 



(1-z) 



l+e 



*(1 " 2) 



l -£-/p-) l+e 
r(2-e)r(-e 



T(2 - 2e) 



[5(£~ - u) - 5{p~ - u)] 

[5{£- - u) - 5(p- - u)] 
(3.8) 



In the first step we used n^Vn = 2 and u n {p)^u n {p) = 2p~ , performed the £ + integral by 
contours and did the £± integral as usual. The £ + integral has the same pole structure as in 
Fig. 5(a)| (except that the double pole at £ + = is now a single pole), which restricts the £~ 
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integral to the finite range < I < p . Expanding Eq. ( |3.S|) in e, using the distribution 
identity in Eq. (|A.3|), we get 



TT 



9(z 
+ S(l-z 



^ 5(1 - z)+ £ (1 - z)z - eCxil - z)z 



1 / 7T 

7 + 1 + £ ( 2 -lT 



2x 1 



6(z)U~ - In j^j [C (l - z)z + 5(1 - z)] - A(l 



+ S(l-z)[2-- 



IT 



z)z 



(3.9) 



where C n (x) = [0(x)(ln n x) /x] + are the usual plus distributions defined in Eq. (|A.2|) . 

In the last step in Eq. (|3.§|), the t~ integral produces an additional 1/e pole in each of 
the two terms corresponding to real and virtual radiation from the two different Wilson line 
contractions. It comes from the singularity at t~ = p~, where the gluon in the loop becomes 
soft. (This soft IR divergence appears as a pole in e because the offshellness only regulates 
the collinear IR divergence here.) The soft IR divergences cancel in the sum of the virtual and 



real contributions, as can be seen explicitly in the first line of Eq. (3J}) where the 1/e poles 
in curly brackets cancel between the two terms. One can already see this in the t~ integral 
in Eq. (|3.8|), because for l" = p~ the two 5 functions cancel so there is no soft divergence in 



the total integral. Thus, in agreement with our discussion in Sec. 2.1, we explicitly see that 



contributions from the soft region drop out in the PDF. As a consequence, the PDF only 
contains a single 1/e pole and correspondingly its RGE will sum single logarithms associated 
with this purely collinear IR divergence. 

Since the gluon in the loop is supposed to be collinear, the soft gluon region must be 
explicitly removed from the collinear loop integral, which is the condition p ^ in Eq. (2.2). 
For continuous loop momenta this is achieved by a zero-bin subtraction. However, since the 
soft region does not contribute to the PDF, it also does not require zero-bin subtractions in 
SCET. (If we were to include separate zero-bin subtractions for the virtual and real contri- 
butions, they would simply cancel each other.) We will see shortly that the situation for the 
beam function is quite different. 



The last diagram with external quarks, Fig. 5(c) 



is 



{qn\Q h q we ^)\q n 



6(l-z)(Zz-l) 



■^5(1-4- 
47r Me 



[X 1 



(3.10) 



Here we used the result for the one-loop on-shell wave-function renormalization with an 
offshellness IR regulator, which is the same in SCET and QCD. 

Adding up the results in Eqs. ([TT]), (^.9[ ), and ( 3.1CQ we obtain for the bare one- loop 
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quark matrix element 



2tt 



In- 



P 3g (^)-£i(l-^)(l + 2 2 



7T 



;i-z)2(i-z) 



where 



P 9<? (z) = £ (l- 2 )(l + ^) + - ( 5(l- Z ) 



0(1 -z 



1 + z 2 



(3.11) 
(3.12) 



is the q — > qg splitting function, see Eq. ( [2.20 ). 

Next, we consider the matrix element of Q q between gluon states \g n ) = \g n (p))- The 



only relevant diagram is shown in Fig. 5(d) 



(gn\QT e {uj)\g n 



4tt 



~g 2 T F 



d d £ (-efa)tT[V»VZ%] 



(2vr) c 



{£ 2 + i0) 2 [(£ - p) 2 + iO] 



5{t 



(3-13) 

Here e = e(p), Vn = Vn (£— p, £) and V£ = V^(£,£ — p). Since the physical polarization vector 
is perpendicular, n ■ e(p) = n ■ e{p) = 0, we only need the perpendicular parts of the collinear 
vertices. The numerator then becomes 



tr 



yuyu/J_ 
' n ' n ^ 



r) 2 (r- P -) = itr 



(p- 



p~ 



-p- 
W + 



t 



\--p- 
1 



l)-\2(n- 



P 



l-z 



(3.14) 



In the last step we used that under the integral we can replace t~ = uj = zp~ and £^_£^_ = 
f^±9jL /(d ~ The remaining loop integral is exactly the same as in Fig. 5(a), so the bare 
one-loop gluon matrix element becomes 



{9n\Qt ve ^)\a n 



(i) _ a s (n)T F nl A nM .A-n/^/^^ 2 



2vr 
2vr 



(z)8(l-z)T(e) 



f 



[1- z)- e (l-2z + 2z 2 -e) 



f)( : ) - In ^ - - *)] - 0(1 - *) }> • ( ■'>■ Lo) 



Here 



P ?9 (z) = 0(1 - z) (1 - 2z + 2z 2 



(3.16) 



is the g ^ qq splitting function from Eq. fl2.20| ) . 

Note that the diagram analogous to Fig. |5(d)| with the two gluons crossed can be obtained 
from Fig. |5(d) by taking p^ — > —p^", which takes z — > — z. The limits resulting from the £ + 
integral are then — 1 < z < or — p~ < uj < 0, and since we require uj > for Q q , this 
diagram does not contribute. The diagram involving the SCET vertex with two collinear 
gluons vanishes because here the £ + integral does not have poles on both sides of the axis. 
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Prom the bare matrix elements in Eqs. and (|3.15| ) we can obtain the renormalization 

of Q q . Taking the quark and gluon matrix elements of Eq. (|2.14| ) and expanding to NLO, 



<gn|Sr e Hkn> 
'da; 
- / u7 

3 



(1) 



E 



^(^^)(^|Qi(^^)kn> ( n^°)(^, / ,)<, n |S i (a;' ) p)| g „>« 



Zl^{z^) + {q n \Q q {oj^)\q n ) {l) 
(i) 



{g n \Q h q are (u)\g. 



(3.17) 



where we used the tree-level matrix elements in Eqs. (2.41) and ( |2.45 ) and Z-j°\z,fi) = 
5ij 5(1 — z). The MS counter terms required to cancel the 1/e poles in the bare PDF matrix 
elements are then 



Zfr) =5(1"*) + 



0(z)P qq (zj, 



1 a s ((x)T F 



2tt 



(3.18) 



Expanding Eq. (|2.15|) to NLO, the one-loop anomalous dimensions are obtained by 

7£(*,a0 = -n-^zl® (j,—a s ((j,) = -2ea s (/i) + /3[a s (/x)], (3.19) 

which with Eq. ( |3. 18 ) yields the anomalous dimension for the quark PDF in Eq. ( |2.19| ). 

= ^^0(z)P qq (z), jf g (z,ri = ^^9(z)P qg (z). (3.20) 

11 7T 7T 

Finally, the renormalized NLO PDF matrix elements, which we will need for the matching 
computation in Sec. 3.3 below, are 



(qn\Qq{^,^)\qn 



(i) _ a a ([i)Cf 



2tt 



{gn\Q q (^,^)\g n ) 



(1) 



g s (fj)T F 
2vr 



(-')<j P gg (z) l n ^+£i(l-z)(l + z 2 ) 

tfCi-^fJ-V) + 0(1-2)2(1-*; 



fl( : ) { P qg (z) [in + ln(l - «)] + 0(1 - 2) 



(3.21) 



3.2 Quark Beam Function with Offshellness Infrared Regulator 

Next, we calculate the bare beam function 5-matrix elements, 



(q n {p)\e(Lj)O h q ^(t,u)\q n (p)) , (g n (p)\e(u)O h q ^(t,u)\g n (p) 



(3.22) 
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Figure 6: One- loop diagrams for the quark beam function. The minus momentum uj is incoming at 
the vertex and the b + momentum is outgoing. Diagram (d) denotes the wave-function contribution. 
Diagrams (b), (c), and (d) have symmetric counterparts which are equal to the ones shown and 
included in the computation. Diagram (f) and the diagram with the gluon connecting both vertices 
vanish. 



to NLO. The corresponding one-loop diagrams are shown in Fig. H. The matrix elements 



are calculated as explained in Sec. |2.5| in Eq. ( 2.5q ): For the virtual diagrams with vacuum 
intermediate state we explicitly insert the vacuum state, while for the real-emission diagrams 
we use Eq. ( 2.571 ). I n the latter case, we first take the Disc, then expand in e to extract the 
UV divergences, and at last take the t' — > limit to isolate the IR divergences into In t' terms. 
Some helpful formulas for calculating the discontinuity and taking the limit t' — > are given 
in App. |A[ 

For the beam function calculation the p + < actually plays a dual role: For the UV 
divergent piece we can treat the calculation as in SCETj, and so p + ~ b + ~ A 2 p~, which 
allows us to explicitly check the structure of the convolution in Eq. (2.22). The renormalized 
result contributes to the matching onto PDFs, matching from SCETj onto SCETn- In the 
matching, — p + <C b + plays the role of the IR regulator, since we are required to use the same 
states as in the PDF calculation. We will see that the IR divergences hit' match up with 
those present in the PDF calculation, and hence drop out in the matching coefficients lij. 

The diagrams in Fig. || have the same Dirac and propagator structure and overall factors 
as the corresponding PDF diagrams in Fig. ^, so we can reuse those parts from the previous 
subsection. The difference compared to the PDF calculation is that for the real-emission 
diagrams, instead of doing the £ + integral by contours, £ + is fixed by the additional 5 function 
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in b + , and since we use time-ordered perturbation theory we must now take the discontinuity. 
This also alters the structure of the remaining £± integral, for which we now use Feynman 
parameters to combine the denominators. After carrying out the £± integration, we will need 
the following two Feynman parameter integrals 



h(A,B,e) 



da [(1 - a) A - aB] 



-l-e 



(-BY 



A- 



da(l -a)[(l -a)A-aB] 



-l-e 



e(A + B) ' 



I l-e 



A- 



(3.23) 



e(l-e)(A + B) 2 e(A + B) 

The first diagram, Fig. 6(a), has real radiation in the final state, so we use Eq. ( 2.57] ) for 
the computation, 



(q n \e(u;)<D h ™ e (t,u;) 



47T 
4-7T 



g 2 c F - 



Disc 



p"(d-2)^ 



(27r) d (^ 2 + i0) 2 [(^-p) 2 + i0] 



<5(r- w)<5(£ + +6 + -p+) 



fl(z)(rf- 2) 
(2vr) 2 z 
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(1-a) 



d d -*lL 

(2vr) d " 2 j$ + (1 - a )A - aB] 3 



a s (n)C F 6(z) 



2ir z 
where we abbreviated 



r(l + e)(e^/i 2 ) e (l - e) 2 [--^Disc t>0 J 2 (A 5, e) 



A = t + t' . 



B 



1 



i. 



,4 + B 



+ f 



(3.24) 



(3.25) 



z z 

Since t' > and z > 0, the only discontinuity in i2(A B, e ) f° r £ > arises from (-B). Using 
Eq. ( |A.7| ) to take the Disc, we obtain 

sinvre 9{B)B l - t 



i i (— B) x ~ e 

--msc t>0 1 2 (A,B,e) = -Disc, >0 e(1 _ e)( ^ + i3)2 



9(t)- 



1 - z 



o{ty 



vre(l-e) (A + 5) 2 ' 
sinvre [(1 - z)t] x ~ 6 z x+t 



'^(l-e) (t + zt') 2 ' (3 ' 26) 
Note that there is only a discontinuity for B > 0, so taking the discontinuity for t > requires 
(1 — z)/z > 0, and since z > we obtain the expected limit z < 1. Since there are no UV 
divergences, we can let e — > 0, and Eq. ( |3.24| ) becomes 



(q n \9(Lu)O b g ^(t,L0)\q n 



(a) _ a s (fx)C F Qf ^ aM ^ 0(i)£ 



0(^)0(1 -z)(l-z)- 



(3.27) 



2vr " v "'" v " y (t + zt') 2 ' 

The above result has a collinear IR singularity for t — > which is regulated by the nonzero 
t' . We can isolate the IR singularity using Eqs. O and flO|) by letting /3 = zt'//u 2 -> 
while holding t = t + zt' fixed 2 , 



lim 



0(t)t 



lim 



t'->0 (i + zt') 2 zt'/n 2 -*) 



6(t-zt') 0(t-zt')zt' 



t 



t 2 



— o^0\ > 



^,/)(ln^ + l) . (3. 



We keep the dependence on £ in our calculation as it will be useful for checking the structure of the 
renormalization in the following subsection. 



- 29 - 



The final result for Fig. 5(a) is thus 



(a) _ a s {n)C f 



2vr 



(z)9(l-z)(l-z) 



{^°(^)-^( ln 



zt 

1? 



+ 1 



(3.29) 

Next, we consider the real-emission diagram in Fig. 6(b) . It corresponds to the 5(£~ —oo) 
term in Eq. Q3.8|). Together with its mirror graph, giving an identical contribution, we obtain 
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(3.30) 



where in the second step we performed the loop integral as before, and in the last step we 
used Eq. ( A. 7 ) to take the discontinuity. As for Fig. 3(a)| , the loop integral produces no 
UV divergence. However, as in the PDF calculation for Fig. |5(b)| , there is a soft gluon IR 
divergence at z — > 1 or £~ — > p~ producing a <5(1 — z)/e IR pole when expanding the last 
factor using Eq. ( |A.3j ). In contrast to the PDF calculation, the soft gluon region must now 
be explicitly excluded from the collinear loop integral. In dimensional regularization with 
an offshellness IR regulator the relevant zero-bin integral is scaleless and vanishes. Thus, 
including the zero-bin subtraction removes the 1/e IR divergence and replaces it by an equal 
1/e UV divergence such that all 1/e poles in the final result are UV divergences. Expanding 
in e, we have 



-_ + ln * ) + £ (l (| 

e p 2 



and taking the same limit as in Eq. (|3~2"8l) to isolate the IR diver gences, 
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(3.32) 
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the final result for Fig. 6(b) 



is 
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(3.33) 
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For the diagram in Fig. 6(c) (and its mirror diagram) we insert the vacuum intermediate 
state between the fields in O q as in Eq. ( p. 58 ) , resulting in a one- loop virtual diagram involving 
a single field. The calculation is exactly the same as for the 5(p~— ui) term in Eq. (3.8) times 
an overall 5(t), 
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(3.34) 



In the last step we expanded in e and took the IR limit. To be consistent we have to use the 
same IR limit in the virtual diagrams as in the real-emission diagrams above, which simply 
turns the overall 6(t) into a 6(t), 



lim 6(t) 

t'-H) 



lim 6(t- zt') = 8(t) 



(3.35) 



As in the PDF calculation, the UV divergence in the loop produces a T(e) and the soft IR 
divergence a T(— e). The latter is converted by the zero-bin subtraction into a UV divergence, 
producing the 1/e 2 pole. The 1/e 2 poles do not cancel anymore between Figs. 3(b) and 
|6(c) as they did for the PDF in Fig. [5(b) , because the phase space of the real emission in 
Fig. 6(b) is now restricted by the measurement of b + via the 6(£ + + b + — p + ). For the same 
reason Fig. 6(a)| has no UV divergence anymore, while Fig. |5(a)| did. The (1/e) lni' terms in 
Eqs. (|3~33|) and @), which are a product of UV and collinear IR divergences, still cancel 
between the real and virtual diagrams, ensuring that the UV renormalization is independent 
of the IR, as should be the case. 

The final one- loop contribution to the quark matrix element, Fig. 6(d) and its mirror 
diagram, comes from wave-function renormalization, 



( qn \9(cu)O h ™ e (t,oj)\ qn ) (d) = 6(t)6(l - z)(Zt 
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(3.36) 

Adding up the results in Eqs. ( 3.29Q , ( 3.33j ), ( 3.34] ) , and ( 3.36| ), we obtain the bare beam 
function quark matrix element at one loop, 
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We now consider the beam function matrix element with external gluons. The corre- 
sponding diagrams are shown in Figs. 5(e) and 6(f) . For Fig. 6(e) , which is analogous to 
Fig. 5(d) , we find 
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(3.38) 



The loop integral and discontinuity are exactly the same as for Fig. 6(a) , The diagram in 
Fig. |6(f) does not contribute to the quark beam function. It can be obtained from Eq. ( 3.38 ) 
by replacing p^ — > —p^, which takes t' —> —t' and z — > —z. Doing so, the only contribution to 
the discontinuity is still from B = —(l + z)t/z for B > 0, which for t > requires —1 < z < 0. 



Hence, Fig. 6(f) does not contribute. Using Eq. fl3.28|) to take t' — v in Eq. ( 3.38] ), we get 
the final result for the bare one-loop gluon matrix element 
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(3.39) 



As for Fig. 3(a), it has no UV divergences because of the measurement of 6 + , which means 



that the renormalization does not mix O q and O g . 



3.3 Renormalization and Matching 



Using the bare matrix elements calculated in the previous section, we can extract the renor- 
malization of O q . We first take t = t + zt' — > t in the bare matrix elements. Then, expanding 
the quark matrix element of Eq. ( [2.22 ) to one-loop order, 
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(3.40) 



= Zf\t^)5{l -z) + {q n \O q (t,u,n)\q n ) {X) , 
we can then read off the MS renormalization constant from Eq. ( |3.37| ) 
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(3.41) 



The fact that the gluon matrix element is UV finite and the UV divergences in the quark 
matrix element are proportional to 5(1 — z) confirms at one loop our general result that 
the renormalization of the beam function does not mix quarks and gluons or change the 
momentum fraction. 
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In Eqs. flO| ) and ( jOll) we used that we already know the structure of the renormal- 
ization from our general arguments in Sec. |2.2| , i.e. that Z q B only depends on the difference 
t — t'. Alternatively, we can also use the dependence on z and the finite dependence on t' via 
t to explicitly check the structure of the renormalization. In this case, we must use the same 
IR limit also for the tree-level result in Eq. ( p. 41 ), which using Eq. (3.35) becomes 



{q n \6{u)O q {t,u,v)\<ln) m = lim S(t) 5(1 - z) = 5(t) 5(1 - z) . 



(3.42) 



Taking Z B (t,t' ,uj/uj' , n) to be a general function of t, t' and uj/lo', we now get for Eq. ( p. 40 ) 
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(3.43) 



In the first step we used Eq. (|1|) and Z q B {0) (t,t' , z) = 5(t - t')5(l - z). From Eq. fl3~37p we 
now find 
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thus explicitly confirming at one loop that Z B (t, t', z, fi) = Z q B (t — t', [i) 5(1 — z). 

The one-loop anomalous dimension for the quark beam function follows from Eq. (|3.41| 
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(3.45) 



It is identical to the one-loop anomalous dimension of the quark jet function. The coefficient 
of Co(t / n 2 ) / ji 2 can be identified as the one-loop expression for — 2T q usp . Thus, Eq. ( |3.45|) 
explicitly confirms the general results in Eqs. ( |2.26| ) and ( |2.35D at one loop. 

Taking the bare matrix elements in Eqs. (|3.37l ) and (|3.39| ) and subtracting the UV di- 
vergences using Eqs. (pUOl) and fljUlp gives the renormalized one-loop beam function matrix 
elements, 
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(3.46) 



For the matching onto the PDFs, we must take t' — > and have therefore set t = t everywhere, 
only keeping t' in the IR divergent hit' terms. 
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Expanding the OPE for the quark beam function, Eq. ( 2.38 ), to one loop, we have 
(q n \6(u)O q (t,u,n)\q n ) (1) 

3 L 
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(3.47) 



Thus, the one-loop matching coefficients, lQ'(t, z, //), are obtained by subtracting the renor- 



malized PDF matrix elements in Eq. ( |3.21| ) from those in Eq. ( |3.46 ). Doing so, we see that 
the hit' IR divergences in Eqs. (|3.21 ) and (|3.46|) precisely cancel, as they must, such that the 
matching coefficients are independent of the IR regulator and only involve large logarithms 
that are minimized at the scale /j, 2 ~ t. The final result for the NLO matching coefficients is 
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4. Numerical Results and Plots 



Including the RGE running in Eq. ( |1.6|) , the full result for the resummed cross section for 
isolated Drell-Yan differential in q 2 , Y, and tb is 
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tg + h 



(4.1) 



where the sum runs over quark flavors ij = {uu, uu, dd, . . .} and the additional contributions 
from the leptonic matrix element are contained in the hard function. Equation ([O]) is valid 
to all orders in perturbation theory. The all-order solutions for the evolution factors in terms 
of the respective anomalous dimensions are given in App. The hard, beam, and soft 
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Table 1: Order counting in fixed-order and resummed perturbation theory. 



functions are each evaluated at their natural scales, — Qi — \Jtb~Q, f^s — t bQi an d 
are then evolved to the common arbitrary scale [i by the evolution kernels Ujj, U % g , and Us, 
respectively. With the one- loop results for the beam function presented above, Eq. (4.1) can 
be evaluated at NNLL order in resummed perturbation theory, which requires the one-loop 
matching corrections, the two-loop standard anomalous dimensions, and the three-loop cusp 
anomalous dimension (see Table |l]). All the necessary ingredients are given in App. If we 
let vb — iO be the Fourier conjugate variable to tb, then the NNLL cross section resums the 
following terms 

In dq2 ^Ydv B ~ ln VB ^ S ln + ^ S ln + as(yOLs ln 
for all integers k > 0. 

In the remainder we will focus on the beam functions, which are the topic of this pa- 
per. Below we compare results for the quark beam function at LO and NLO in fixed-order 
perturbation theory as well as at NLL and NNLL in resummed perturbation theory. Our 
conventions for the a s loop counting are given in Table |]. To evaluate the required convolu- 
tions of plus distributions at NNLL we use the identities from App. B of Ref. [27]. We always 
use the MSTW2008 @ parton distributions at NLO for a s (m z ) = 0.117 and with two-loop, 
five-flavor running for a s (fj,). The uncertainty bands in the plots show the perturbative un- 
certainties, which are estimated by varying the appropriate scales as explained in each case. 
They do not include the additional uncertainties from the PDFs and a s (mz)- 

The order of the running of a s (fj,) deserves some comment. Working consistently to NLO 
in the matching corrections requires us to use NLO PDFs, for which the two- loop running 
of a s was used in Ref. [p2^ . On the other hand, the double-logarithmic running of the hard 
function and beam functions at NNLL requires the three- loop running of a s , which poses a 
slight dilemma. Ideally, we would need NLO PDFs using three- loop running for a s (fi), which 
as far as we know is not available. The numerical difference between a s run at two and three 
loops is very small, at most 2%. Hence, we use the following compromise. To be consistent 
with our PDF set, we use the above a s (mz) and two-loop, five-flavor running to obtain the 
numerical value of a s at some required scale, and to be consistent with the RGE, we use the 
two- and three-loop expression for the QCD (3 function in the RGE solutions at NLL and 
NNLL. (For simplicity we use the same NLO PDFs and a s also at NLL.) 

To illustrate the importance of the various contributions to the quark beam function, 
we also consider the beam function in the threshold limit and without the gluon mixing 
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Figure 7: The u (left) and u (right) beam functions at the hard scale = x 7 TeV at LO, NLO, NLL, 
and NNLL, integrated up to i max = (xe~ 2 7TeV) 2 . The bands show the perturbative uncertainties 
estimated by varying fin for the fixed-order results and the matching scale [x 2 B ~ i m ax for the resummed 
results, as explained in the text. 



contribution. In the threshold limit we only keep the terms in Eq. (3.48) which are singular 
as z — > 1, 
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(4.3) 



The gluon mixing term I qg contains no threshold term (which reflects the fact that in thresh- 
old Drell-Yan the gluon PDF does not contribute). For the result without the gluon mixing 
contribution we keep the full X qq but set I qg to zero, which corresponds to adding the remain- 
ing non-threshold terms in X qq to the threshold result. In the plots below, the results in the 
threshold limit are shown by a dotted line and are labeled a x — > 1" , and the results without 
the gluon contribution are shown by a dashed line and are labeled "no g" . The full result, 
including both X qq and I qg , is shown by a solid line. Hence, the size of the non-threshold 
terms in I qq , and therefore the applicability of the threshold limit, is seen by the shift from 
the dotted to the dashed line, and the effect of the gluon mixing is given by the shift from 
the dashed to the solid line. 

To be able to plot the beam function as a function of the momentum fraction x including 
the virtual terms proportional to S(t), we consider the integral over t up to some maximum 



Bi(t max , x, /i) — / dt Bi(t, x, jJL)6(t 



t) 



(4.4) 



where Bi(t,x,fi) is given by Eqs. ( 2.38j ) and ( 2.28j ), In the plots, we always choose t 



(xe" 2 7TeV) 2 , which one should think of as t r 



(e-y cut xE c 



Hence, this choice of 
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Figure 8: The u (left column) and d (right column) beam functions at the beam scale y? B ~ t max at 
LO and NLO, integrated up to i max = (xe~ 2 7TeV) 2 . The top row shows the functions times x. The 
bottom row shows the relative differences compared to the LO result. Also shown are the NLO beam 
functions in the threshold limit (dotted) and without the gluon contribution (dashed). The bands 
show the perturbative scale uncertainties as explained in the text. 



^max corresponds to a rapidity cut y cut = 2 for E cm = 7TeV or equivalently y cut = 2.4 for 
E C m = lOTeV. This is motivated by the upper bound y cut = y B ut ± Y, which follows from 
the factorization theorem Eq. (4.1) when we integrate tb < exp(— 2y c B at ). 

Figure ^ shows the integrated u and u-quark beam function xBi(t max , x, hh) evaluated 
at the hard scale [in = Q = x7TeV. For the fixed-order results at LO (lowest gray band) 
and NLO (wide green band), the bands are obtained by varying hh by factors of two, since 
this is the scale at which the perturbation series for the matching coefficients in Eq. ( |2.38| ) is 
evaluated. At LO, the resulting variation is entirely due to the scale dependence of the PDF. 
At NLO the sizeable variation indicates the presence of large single and double logarithms of 
tmax/Q 2 when the fixed order beam functions are evaluated at \ih- 

For the resummed results at NLL (blue band) and NNLL (orange band) the beam function 
OPE, Eq. ( p.38| ), is evaluated at the beam scale fi B ~ i m ax> and the beam function is then 
evolved to fin using its RGE, Eq. (|J|). In this way, the large logarithms of ^ B /^? H ~ 
^max/ M# are resummed. Here the bands correspond to perturbative uncertainties evaluated 
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Figure 9: The u (left column) and d (right column) beam functions at the beam scale. The meaning 
of the curves is analogous to Fig. M. 



by varying the matching scale \ib while keeping fixed. The dependence on the scale \xb 
cancels between the fixed order perturbative result for the beam function and its evolution 
factor, up to higher order corrections. An estimate for these higher order corrections is given 
by the NLL and NNLL bands. These uncertainty bands show the minimum and maximum 
variation in the interval \/t max /2 < \ib < 2y/t ma _ x (which due to the double-logarithmic series 
do not occur at the edges of the interval) with the central value given by the center of the 
bands. The NLL result is close to the NLO result, showing that the large logarithms make up 
by far the biggest part in the NLO corrections. Consequently, the corrections from NLL to 
NNLL are of reasonable size and within the NLL uncertainties. Hence, for the beam function 
at the hard scale, fixed-order perturbation theory is not applicable. Resummed perturbation 
theory is well-behaved and should be used. 

To study the perturbative corrections to the beam functions in more detail, we consider 
them at the scale y? B ~ t max , where there are no large logarithms and we can use fixed-order 
perturbation theory. The u and d beam functions at LO and NLO are shown in Fig. ||, and the 
u and d beam functions in Fig. ^. The top rows show xBi(t m3 ^j x, hb)- The bottom rows show 
the same results but as relative corrections with respect to the LO results. At LO, the only 
scale variation comes from the PDFs and the minimum and maximum variations are obtained 
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for hb = {\/*max/2, with the central value at hb = v^max- For the NLO results, 

the maximum variation in the range y/t max /2 < ^b < 2i/imax is approximately attained for 

= {0.7-v/tmax, 2.0^t max } and the corresponding central value for /ig = 1.4\/f max . To be 
consistent we use the same central value /is = 1.4\/i max for the NLO results in the threshold 
limit and without gluon contribution. In all cases the NLO perturbative corrections are of 
0(10%) and exhibit reasonable uncertainties. 

The integration limits x < £ < 1 in the beam function OPE, Eq. ( |2.38| ), force z = x/£ — > 1 
in the limit x — > 1. Hence, the threshold terms in Eq. (4.3) are expected to dominate over the 
non-threshold terms at large values of x. This can be seen in Figs. |8]and[9|, where the threshold 
results shown by the dotted lines approach the full results towards large x values where the 
beam functions vanish. For the quark beam functions in Fig. |8| away from the endpoint, 
x < 0.5, the threshold corrections give a poor approximation to the full NLO corrections. For 
the antiquark beam functions in Fig. |9[ the threshold result turns out to be relatively close to 
the full result even for small x. However, the reason for this is a relatively strong cancellation 
between the non-threshold terms in the quark and gluon contributions X qq and X qg at one 
loop. As shown by the result without the gluon contribution (dashed lines) the non-threshold 
terms in the quark and gluon contributions each by themselves are of the same size or larger 
than the threshold contributions. Note also that for the d beam function the threshold result 
approaches the no- gluon result rather than the full result at large x. A similar but less strong 
cancellation can also be observed at small x in the quark beam functions. These appear to 
be accidental cancellations, which depend on both the relative size of the (anti)quark and 
gluon PDFs as well as the relative size of the non-threshold terms in I qq and I qg . Therefore 
one must be careful when applying the numerical dominance of the threshold terms to cases 
where it is not explicitly tested. 

It has been argued |33], 34] that the steep fall-off of the PDFs causes a systematic en- 
hancement of the partonic threshold region z — > 1 even away from the hadronic threshold 
limit x — > 1. This likely explains why the threshold terms in Figs. |8| and ^ start to dominate 
already close to the x values where the PDFs are close to zero, rather than strictly near 
x = 1 |J5|. However, our results show that the same arguments do not apply in the relevant 
region of x where the PDFs and beam functions are substantially nonzero. 



5. Conclusions 



At the LHC or Tevatron, the appropriate description of the initial state of the collision depends 
on the measurements made on the hadronic final state. The majority of measurements trying 
to identify a specific hard interaction process do so by finding a certain number of central jets, 
leptons, or photons that are distinguished from energetic initial-state radiation in the forward 
direction. These measurements effectively probe the proton at an intermediate beam scale 
Hb *C Q and the initial state is described by universal beam functions. The beam functions 
encode initial-state effects including both PDF effects as well as initial-state radiation forming 
an incoming jet around the incoming hard parton above [ib- 
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We have discussed in detail the field-theoretic treatment of beam functions using SCET. 
We discussed their renormalization properties and showed that they satisfy an RGE with the 
same anomalous dimension as the jet function to all orders in perturbation theory. The beam 
function RGE determines the evolution of the initial state above [Xb- It resums a double 
logarithmic series associated to the virtuality t of the incoming parton, while leaving the 
parton's identity and momentum fraction x unchanged. 

We gave a general discussion of the operator product expansion for the beam functions 
that allows us to match them onto PDFs fj(£, lib] convoluted with matching coefficients 
Iij(t, x/£, Mb). The latter encode the effects of the initial-state radiation and are perturba- 
tively calculable at the scale lib- We performed this matching at one loop for the quark beam 
function onto quark and gluon PDFs. Our calculation explicitly confirms at one loop that 
the quark beam function contains the same IR singularities as the PDFs, and this required a 
proper handling of zero-bin subtractions. 

In Sec. ||, we presented an explicit expression for the resummed beam thrust cross section 
for Drell-Yan production, pp — > X£ + £~, with the necessary ingredients for its evaluation at 
NNLL collected in App. IDL An analysis of the cross section at this order is left to a separate 



publication [36]. Here, we discussed in detail numerical results for the quark beam function 
at NLO and NNLL. The gluon beam function is important for Higgs production at the LHC. 
The one-loop matching of the gluon beam function onto gluon and quark PDFs is discussed 
in a separate publication and used to calculate the Higgs production cross section for small 



beam thrust at NNLL [37]. Another application is to define a px dependent beam function 
to study the pr-spectrum of the Higgs [Q . 

So far, effects of strong initial-state radiation that go beyond the inclusive treatment via 
PDFs have only been studied using Monte Carlo methods and models for initial-state parton 
showers. The physical picture behind the beam function and initial-state parton showers are 
in fact in close correspondence. Beam functions and the beam thrust factorization theorem 
provide a complementary field-theoretic approach to study these effects analytically. Hence, 
they provide a crucial tool to obtain an accurate description of the initial state, which is 
mandatory to obtain precise and realistic theory predictions for the LHC. 
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A. Plus Distributions and Discontinuities 

The standard plus distribution for some function g(x) can be defined as 

[6(x)g(x)] + = lim A [0( x - p) G{x)} with G(x) = J dx' g{x') , (A.l) 
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satisfying the boundary condition dx [9(x)g(x)]+ = 0. Two special cases we need are 

~9(x)ln n x~ 



C n (x) 
C r '{x) 



lim 

P^yO 



x n + 1 



In addition, we need the identity 

0(s) 



lim 

/3->0 



1 



-^ l + S(x-p) 



X 



n 



(A.2) 



»l+£ 



-5(x)+£ (x)-e£i(x) + O(e"), 



the Fourier transform 



U{x) = - f^-e^ ln[i(y-i0)e^], 



(A.3) 
(A.4) 



and the two limits 



lim 

p-to 



x 2 



7T _ 



b 



/3^0 X 2 



(A.5) 



Away from x = these relations are straightforward, while the behavior at x = is obtained 
by taking the integral of both sides. General relations for the rescaling and convolutions of 
C n (x) and C^{x) can be found in App. B of Ref. [27]. 
The discontinuity of a function g(x) is defined as 



T>\sc x g(x) = lim [g{x + i/3) - g(x - i/3)l 

P-tQ 



(A.6) 



If we are only interested in the discontinuity in some interval in x, we simply multiply the right- 
hand side with the appropriate 9 functions, as in Eq. ( p. 55 ). If g{x) is real then ~D\sc x g{x) = 
2ilm g(x + iO). Two useful identities are 



i _ 1 
— Disc 
2vr 



ra+1 



nl 



l 

2^ 



Disc x (-x) n - £ = (-l) n - lSm7re 



7T 



9(x)x n - e . (A.7) 



To derive the last identity, note that (— x— i0) n e = exp[(n— e) ln(— x— iO)] = \x\ n e exp[— m(n- 
e)6(x)], so taking the imaginary part gives Im(— x — i0) n_e = (— l) n sin(7re) 9{x) x n ~ e . 



B. Renormalization of the Beam Function 



In this appendix we derive the general structure of the beam function RGE in Eq. ( p. 25 ) 
to all orders in perturbation theory. The two essential ingredients will be the known all- 
order renormalization properties of lightlike Wilson lines [39, 4(| 41, 42] and the factorization 
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theorem for the isolated pp — > XL cross section, where X is the hadronic and L the non- 
hadronic final state. In Ref. Q we proved that to all orders in perturbation theory and leading 
order in the power counting this cross section factorizes as 



da 



d^Wd^d-B^ 

-.2 



Hij{q 2 ,Y, fi) f dk+ dk+ Sf[ emi (k+, fc+ /*) 
x q 2 Bi[uj a (B+ - k+),x a ,fi}B j [u b {B+ - k+),x b ,n] . (B.l) 



The sum over ij runs over parton species ij = {gg, uu, uu, dd, dd, . . .}. The soft function does 
not depend on the quark flavor, and its superscript only refers to the color representation. 
The variables q 2 and Y are the total invariant mass and rapidity of the non-hadronic system 
L, x a ,b = \fq 2 e ±Y / E cra and uj a> b = x a ^E cra . The hadronic variables B^ b are the hemisphere 
plus momenta of the hadronic final state X with respect to the directions n a and rib of the 
incoming protons. Their precise definition will not be relevant for our discussion. 

The three ingredients in Eq. flB.ip are the renormalized hard, beam, and soft functions, 
Hij(q 2 ,Y, fj,), Bi(t,x,fi), -(k£ , k^ , fi) . Their dependence on the renormalization scale \i 
must cancel in Eq. ( |B.1[) , because the cross section must be [i independent. The structure of 
the RGE for the hard and soft functions thus uniquely determines the allowed structure of 
the beam function RGE. 

The hard function is a contraction between the relevant leptonic matrix element squared 
and the square of the Wilson coefficients of the color-singlet qq and gg local SCET currents 

o# = xt,-. a , o% = v^c.-^a ®Z-« b ± ■ ( B - 2 ) 

where a and (3 are spin indices. In each collinear sector, total label momentum and fermion 
number for each quark flavor are conserved. Thus, the currents cannot mix with each other 
and are multiplicatively renormalized. Furthermore, RPI-III invariance implies that the RGE 
for the currents can only depend on q 2 = uj a uJb- The renormalization of these SCET currents 
also does not depend on their spin structure, so the RGE for the hard function must have the 
same structure as for the currents. Therefore, to all orders in perturbation theory we have 
(with no sum on ij) 

^H lJ (q 2 ,Y,ri)= 1 %(q 2 , t i)H ij (q 2 ,Y,f,). (B.3) 

Next, the incoming hemisphere soft function, emi (fc^") k^, /i), is given by the vacuum 
matrix element of incoming soft lightlike Wilson lines along the n a and n& directions. In 
position space, 

SLjUaiVb ,M) = / dktdk^e-^^y 2 Sii cmi (kt,K^) (B.4) 
has two cusps, one at spacetime position and one at y = y^n a /2 + y^rib/2. The renormal- 



ization properties of lightlike Wilson lines with cusps [39, 40, 41, 42 1 then imply that to all 
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orders in perturbation theory, 



(B.5) 



Js(y a ,y b »m) = 2r cus P (a s ) 



In i 



iO 



-/xe 



In i 



iO 



-/xe 



+ 7^) 



where is the cusp anomalous dimension for quarks/antiquarks or gluons, and 7^ [a s (/z)] 

and r* usp 



a s (fj,)] depend only indirectly on fi via a s ([i). Dimensional analysis and RPI-III 
invariance imply that the single logarithm multiplying 2r* usp scales like ln(y~y^~// 2 ). (The 
additional dimensionless factors are chosen for convenience. Any change in them can be 
absorbed into 7^ (a s ).) The correct overall sign and iO prescription for the logarithms can be 
deduced from the explicitly known one- loop result [43], • 

Taking the Fourier transform of the cross section in Eq. (BJ) with respect to and 
and differentiating the result with respect to \i yields 



= n 
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2uj a 



x a ,n)B j (^^-,x b ,n 



^ihcmiVya 
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B, 



( Va 
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,x a ,fj,)Bj\^-,x b ,iJ, 



(B.6) 



The factorization theorem for the cross section neither depends on the choice of L, which 
affects the form of Hij for different ij, nor the type of the colliding hadrons. This implies 
that each term in the sum over ij must vanish separately. (For example, choosing Drell-Yan, 
L = £ + £~, there is no contribution from ij = gg, so the quark and gluon contributions are 
separately zero. Then, by assigning arbitrary electroweak quark charges, the contribution 
from each quark flavor must vanish separately. Finally, the ij = qq and ij = qq contributions 
for a single quark flavor q must vanish separately by choosing various different incoming 
hadrons.) Therefore, the RGE for the product of the two beam functions is 

d 



7^(w a w 6 ,^) + y(y, 



1 ' ' (a, a > y b 



d/i 



Bi 



0. 



(B.7) 



which shows that the beam functions in position space renormalize multiplicatively and in- 
dependently of x a;b . The RGE for each individual beam function can only depend on the 
RPI-III invariant y~ /2uj and obviously cannot depend on the variables of the other beam 
function. Hence, we find that to all orders in perturbation theory 



X,fJ, 



(B.8) 



which is the result we set out to prove in this Appendix. Using Eq. ( ]B.8j ) together with 
Eq. ([B.7D, the anomalous dimensions must satisfy the consistency condition 



= jfiiuaut,, n) + 7^ (y, 



+ 1b 



(Vh_ 

\2uj b 



(B.J 
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Given the form of 7^ in Eq. (B.E), it follows that the anomalous dimensions are given to all 
orders by 



IB 



\2u' 



2Y 
2T 



, \ , U a UJb ij , . 
h 1 



cuspK) Mi 
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iO 



2w 



7s = -7^(a s ) - 7s (a s ) - 7s 



(B.10) 



Taking the Fourier transform using Eq. (| 
become 



2r 



cusp 



, the momentum-space anomalous dimensions 



Uj^-)5{kt) + 5{ki) l -C (^ 
fj, \ fx / fi \ n 



7 i(t,M) = - 2r cusp(« S )^ ^J +7^K)5(t). (B.ll 
The same all-order structure of the soft anomalous dimension as in Eq. ( |B.11| ) was ob 



tained in Ref. [26] for the hemisphere soft function with outgoing Wilson lines in e + e — > 2 
jets using analogous consistency conditions. In fact, the hard SCET currents here and there 



are the same and in Sec. |2.2| we proved that the anomalous dimensions for the beam and jet 
function are the same, 7^ = jj. Hence, the hemisphere soft functions with incoming and 
outgoing Wilson lines have in fact identical anomalous dimensions to all orders. 

C. Matching Calculation in Pure Dimensional Regularization 

Here we repeat the NLO SCETj to SCETn matching calculation from Sec. [3] using dimensional 
regularization for both the UV and IR. Since we only change the IR regulator, the final results 
for the matching coefficients lij(t,z,fi) should not be affected. 

In pure dimensional regularization all the loop diagrams contributing to the bare matrix 
elements of Q q vanish, since by dimensional analysis there is no Lorentz invariant quantity 
they can depend on. Hence, including the counter terms in Eq. ( p. 18 ) to subtract the UV 
divergences, the renormalized matrix elements consist of pure IR divergences with opposite 
signs to the UV divergences, 



(q n \Q q (u,[i)\q ri 

(9n\Qq(v,V)\9n 



1 a s (n)C F 

e 2ir 
1 a s (fJ-)T F 
e 2vr 



0(z)P qg (z). 



(C.l) 



This shows explicitly that the conventional MS definition of the PDFs in QCD, which also 
yields Eq. ( |C. 1| ) , is indeed identical to the SCET definition used in our OPE for the beam 
function. 



Considering the beam function matrix elements, the bare results for Figs. 6(c) and 6(d) 



now vanish, because their loop integrals are again scaleless. For the remaining diagrams we 
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can reuse the intermediate results from Sec. 3.2 before carrying out the Feynman parameter 
integrals and taking the discontinuity. Setting t' = the denominator in the Feynman 
parameter integrals in Eq. ( 3.25| ) becomes (1 — a)A — aB = t(l — a/z). In this case it easier 
to carry out the integral after taking the discontinuity. The discontinuity we need is 



2vr 



-Disc 



*>o 
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" 1_(E _ sin7re 9(t) /a Af a \" ~ ~' 



(C.2) 



where we used Eq. ( |A.7| ). Since we require z > 0, the first 9 function becomes 9 (a — z), and 
so we have 



i . s sin7re 9(t) f 1 , , /a \ -l-e 
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For Fig. 3(a) , using Eqs. ( 3.24 ) and ( |C.3j ) we obtain 



9(1 - z)z 1+t (l - z)- e , 

1- z)z 1+e (l- z) 1 - 6 
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where in the last step we used Eq. (|A.3|) to expand in e. For Fig. 6(b), we start from the 



third line in Eq. (|3.30|) and using Eqs. ( |C.3|) and (|A.3|) we get 
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Adding up Eqs. ( p.4|) and ( |C.5| ), the bare quark matrix element in pure dimensional regular- 
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ization becomes 
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We can now proceed in two ways to obtain the matching coefficient I qq (t, z, /j,). 

First, we can subtract 8(t) times Eq. (|C.1[) from Eq. ( C.6 ) to obtain the bare matching 
coefficient. This simply removes the (l/e)5(t)P qq (z) in the first line of Eq. flC.6|) . Assuming 
that the IR divergences between the PDF and beam function cancel (and including the 
vanishing zero-bin) the remaining poles in the first line are of UV origin and determine the 



necessary MS counter term, reproducing our previous result for Z q B (t,[i) in Eq. ( 3.41 ). 

Alternatively, we can use our general result that the beam function has the same renor- 
malization as the jet function. In this case, we subtract the one-loop counter term for O q are in 



Eq. ( 3.41| ) (which is already known from the jet function's renormalization) from Eq. ( |C.6| ) to 
obtain the renormalized quark matrix element, which equals Eq. ( p.6| ) without the [...]<5(1 — z) 
term in the first line. The remaining 1/e pole must then be of IR origin, so we again have 
an explicit check that the IR divergences in the beam function match those of the PDF in 
Eq. ( p.l| ). Either way, the finite terms in the last two lines of Eq. (C.6) determine the renor- 
malized matching coefficient I qq (t, z, n), which agrees with our previous result in Eq. fl3~48[) . 

For the gluon matrix element, Fig. |6(f) again does not contribute. For Fig. 6(e)| , starting 
from the third line of Eq. ( |3.38| ) , we find 
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The same discussion as for the quark matrix element above can be repeated for the gluon 
matrix element. The (1 / e)5(t)P qg (z) term matches the IR divergence in the PDF in Eq. ( |C . 1| ) . 
Since there are no further poles, no UV renormalization is required and the quark and gluon 
operators do not mix. The finite terms in Eq. (|C.7[) then determine the matching coefficient 
T qg (t, z, /i), reproducing our previous result in Eq. fl3,48|) . 
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D. Perturbative Results 



In this appendix we collect perturbative results relevant for the Drell-Yan beam thrust cross 
section in Eq. fl4.1|). 



D.l Fixed-Order Results 

The one-loop Wilson coefficient from matching the quark current from QCD onto SCET was 
computed in Refs. |44, |45| ], 
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in agreement with the one- loop quark form factors. The hard function is given by the square 
of the Wilson coefficient Pi 



(l-m|/g 2 ) 2 + m|r|/^ J 1 W ' ^ 1 ' 

(D-2) 

where we included the prefactor from the leptonic matrix element, Q q is the quark charge 
in units of |e|, v^ q and a^ q are the standard vector and axial couplings of the leptons and 
quarks, and mz and Y z are the mass and width of the Z boson. 

As discussed in Ref. ]!]], the one-loop result for the beam thrust soft function can be 
extracted from the one-loop incoming hemisphere soft function [43, [2(J, yielding 

a s {\i) C F 



S B (k + ,n) = 6(k + ) + 



2tt 
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(D.3) 



Our one-loop results for the matching coefficients in the beam function OPE in Eq. (2.38) 
are given in Eq. ( 3.48] ). 

D.2 Renormalization Group Evolution 

The RGE and anomalous dimension for the hard Wilson coefficients are [44, 45 1 



The anomalous dimension for the qq hard function in Eqs. flB.3j ) and ( B.lOj ) is given by 
Jh (q 2 > L 1 ) = 2Re[7^((7 2 , fi)]. The expansion coefficients of rc U sp(«s) and 7^(a s ) are given 
below in Eqs. ( D.15| ) and (D.17). The solution of the RGE in Eq. (D.4) yields for the 
evolution of the hard function 
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where the functions K^(fio, fx), r/p(^o, m) and K 7 are given below in Eq. ( D.ll ). 
The beam function RGE is [see Eqs. and ( ^26|) ] 

d f 

IJ--^Bi(t,x,n) = J dt'y B (t-t',n)Bi(t',x,n) , 

7b(*,m) = - 2r cus P («s) ^0 (^) +7bKH(0 , 
and its solution is 0, g|, [§, |7| [see Eq. ( pT28| )l 

Bi(t,x,fi) = / dt' Bi(t - t',x,Mo) U l B (t',no,fi) , 
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The beam thrust soft function is given in terms of S^emi by 
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Its RGE is easily obtained by integrating Eqs. ( |B.5| ) and ( p. 11 
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7s(^ + , /x) = 4r« usp (a s ) -£ (— ) + 7s(a s ) , 7 s(« s ) = -2 7 ^(a s ) - 2 7 « (a s ) 

whose solution is completely analogous to Eq. ( p. 7] ), 
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The functions K^(/io,fj,), rj^(fiQ, m), K 7 (mojM) m the above RGE solutions are defined as 
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Expanding the beta function and anomalous dimensions in powers of a s , 
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their explicit expressions at NNLL are (suppressing the superscript i on K^, r]\ and T l n ), 
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Here, r = a s (fj,) / a s (fJ>o) and the running coupling is given by the three-loop expression 
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where X = 1 + a s (//o)pb l n (A t /A t o)/(2vr). As discussed in Sec. in our numerical analysis we 
use the full NNLL expressions in Eq. ( p,13 ), but to be consistent with the NLO PDFs we 
only use the two- loop expression to obtain numerical values for a s (n), hence dropping the 
(3 2 and /3 2 terms in Eq. flD.14Q . (The numerical difference between using the two-loop and 
three- loop a s is numerically very small and well within our theory uncertainties.) Up to three 
loops, the coefficients of the beta function and cusp anomalous dimension [40, ^] in 

MS are 
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The MS anomalous dimension for the hard function can be obtained [48, 49] from the IR 
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divergences of the on-shell massless quark form factor which are known to three loops |3(| , 
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Denoting ■y'j the coefficient of the 5(1 — z) in the quark PDF anomalous dimension, Eq. ( [2.17D 
(which gives the non-cusp part of the anomalous dimension in the threshold limit z — > 1), 
the factorization theorem for DIS at threshold implies that 2"f q H (a s ) + 7j(a s ) + ji(a s ) = 0, 
which was used in Ref. [49 to obtain 7j at three loops from the known three- loop result for 



[29]. As we showed in Sec. 2.2, the anomalous dimension for the beam function equals 
that of the jet function, 7# = 7j, so the three-loop result for 7 J together with Eq. (DAT) 
yields the non-cusp three-loop anomalous dimension for the beam function, 



7b 



7b l 



7b 2 



6C F , 
C F 
2C F 



/146 

It " oul > 3 

1 52019 841vr 2 82vr 4 



8OC3 )C A + (3 - 4vr 2 + 48Cs)Cf + 



/121 2vr 2 
f + 



2056C 3 



V 162 
/ 151 



+ 



+ 



+ 



+ 



81 
205vr 2 



27 

247vr 4 



V 4 
29 

T 



9 



+ 3vr 2 + 



8^ 4 



7739 325 



54 
3457 
324 



+ 



+ 



81 

57T 2 

9 



vr 2 + 



135 
+ 68C3 - 

6177T 4 



+ 
y 

+ 844C3 + 



88tt 2 C 



V 9 

3 



9 

8tt 2 C 3 



3 

+ 232C 5 



270 



+ 



I6C3 



3 

16^Cs _ 
3 

_ 1276(3 
9 

/ 1166 
V 27 



+ 120(5 )C A C F 



240Cs )C F 



Svr 2 
9 



417T 4 
T35 



+ 



(D.18) 



At NNLL, we only need the one- and two- loop coefficients of 7^ and 7^. The three-loop coef- 
ficients, 7^2 and 7 B2 ' are gi ven here for completeness. They are required for the resummation 
at N 3 LL, where one would also need the four-loop beta function and cusp anomalous dimen- 
sion, the latter of which is has not been calculated so far. In addition, the full N 3 LL would 
also require the two-loop fixed-order corrections, which are known for the hard function, but 
not yet for the beam and soft functions. 
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